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INTRODUCTION 

The principal goal of this project was to develop and utilize successively 
improved versions of a global 3-dimensional dynamical-chemical model of the 
stratospheric ozone layer# The major accomplishments are described in the 
following papers: 

A three-dimensional dynamical-chemical model of atmospheric ozone. _J. 
Atmos. Sci ., 32 , 170-194, (1975) (D. Cunnold, F. Alyea, N. Phillips, R. 
Prinn) . 

Stratospheric ozone destruction by aircraft-induced N0 X * Science , 188 , 
117-121, (1975) (F. Alyea, D. Cunnold, R. Prinn). 

Stratospheric distributions of odd nitrogen and odd hydrogen in a two- 
dimensional model. J. Geophys. Res ., 80 , 4997-5004, (1975) (R. Prinn, F. 
Alyea, D. Cunnold). 

The impact of stratospheric variability on measurement programs for minor 
constituents. Bull. Amer. Met. Soc ., 57 , 686-699, (1976) (R. Prinn, F. 
Alyea, D. Cunnold). 

The dependence of ozone depletion on the latitude and altitude of injection 
of nitrogen oxides by supersonic aircraft. Amer. Inst. Aero, and Astro. 
Journal , 15 , 337-345, (1977) (D. Cunnold, F. Alyea, R. Prinn). 

On the radiative damping of atmospheric waves. J. Atmos. Sci., 34, 1386- 
1401, (1977) (R. Prinn). 

Photochemistry and dynamics of the ozone layer. Annual Rev. Earth Plan . 

Sci . , 6_, 143-174, (1978) (R. Prinn, F. Alyea, D. Cunnold) . 

Preliminary calculations concerning the maintenance of the zonal mean ozone 
distribution in the Northern Hemisphere. Pure and Applied Geophys ., 118 , 
329-354, (1980) (D. Cunnold, F. Alyea, R. Prinn). 
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The above papers described a "6-wave" version of our 3D model. Over the 
past several years we have also been developing a very much improved " 18-wave" 
version. The major differences between the 18-wave and 6-wave models are a 
follows : 

(a) Uses a horizontal diffusion coefficient, A^, to close off the strato- 
spheric jets. 

(b) Non-zonal tropospheric forcing in the tropics (i.e., the annual or even 
components) has been neglected. Zonal heating terms have been reevaluated 
and modified slightly in the troposphere, primarily in the P°(p) component. 


(c) The h profile has been modified in the troposphere to be in closer agree- 
ment with Trenberth. This has the effect of putting the non-zonal forcing 
at slightly higher altitudes and of increasing the rate of dissipation of 
available potential energy in the troposphere relative to some 6 wave compu- 
tations (but not run 12). 

(d) The coefficient of surface friction has been increased from 1.6 to 2. 

(e) The static stability has been returned to its ’'original" value in the tropo- 
sphere. This profile was not used in any reported 6 wave run and is more 
stable than that used in either run 12 or 17. The global mean temperature 
has also been increased by a few degrees at levels 24, 25 and 26. 

(f) The vertical diffusion coefficient is set to 100 cm 2 /sec throughout the 
stratosphere. 

(g) The non-linear Jacobian terms in the model are now evaluated using transform 
methods rather than the original interaction coefficient technique. 


This 18-wave model has now been completed and is described fully in the 
following pages. 
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] . Basic dynamical equations and coordinate system. 


The horizontal coordinate system will be longitude (positive eastward) and 
latitude, denoted by X and <j>. This dependence will be represented in spherical 
surface harmonics, except that certain terms, such as part of the heating and 
photochemistry will be evaluated point-wise at selected values of X and <f>. In 
the vertical direction pressure (p) will be used as a coordinate with finite- 
differences being employed. These pressure levels will be distributed at equal 
intervals of log P in order to give roughly equal intervals in height. We 
define 

P = p 4- (100 cbar) 

Z = -InP, P = e" Z 

From the hydrostatic relation dp = -pgdz and p = p/RT, we have 

dZ = -f E -|r dz (1 - 2) 

The vertical levels will be separated by a uniform value of AZ. To the extent 
that the temperature T is approximately uniform at near surface values, a change 
of one in Z corresponds to a height change of the order of 7 km. The bottom of 
the atmosphere will, for simplicity, be taken at Z = 0, i.e. , at p = 100 cb in- 
stead of at the conventional sea-level pressure of 101.325 cb. The top of the 
"atmosphere" will be artificially set at Z = Z^. corresponding to a geometric 
height of about 70 km. 

The dynamical system not only assumes hydrostatic balance, but also a 
"quasi -geostrophic balance" in the horizontal equations of motion. Because we 
must consider global processes over the entire sphere, this balance must allow 
for complete variability of the coriolis parameter f: 
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f = 2ft sincf) 

ft = 7.292 x ID’ 5 rad sec 


The quasi-geostrophic balance in question is obtained as follows (Lorenz, 

Tel lus , 1960, P. 364). First, we divide the horizontal velocity v into a non- 

A 

divergent part k x Vi() given by a stream function ip and a divergent part 
given by a velocity potential x : 

v = k x 7i|i - Vx (1.4) 


If the eastward and northward components of v are represented by u and v and a 
is the radius of the earth, this is equivalent to 


u - 


fdA 


a COSfcfr- - - 


cM 

34 > 


J lx. 


a cos4> 3A 


v = a 


d<t> 


1 


dt a cos4> 3 A a 3<j> 


li _ 1 

IX a 


9X 


(1.5) 


The vertical component of relative vorticity, ?, and the horizontal divergence 
of "v are related to ip and x by 


C = k • curl v = V z ip; div v = - V 2 x (1.6) 

where V 2 is the horizontal Laplacian operator on the sphere. 

The condition of the quasi-geostrophic balance is 

V • fV\p = gV 2 z (1.7) 

where g is gravity and z is the height of a constant pressure surface. [Unless 
noted otherwise, all partial derivatives with respect to X, -+> , and t (time) are 
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r 


carried out at constant pressure (or Z)]. The hydrostatic relation. 


3z _ 1 _ RT 

9 3p " " p P 

or 


g 


32 r 
3Z " 


RT 


enables (1.7) to be rewritten as 


V 


fvll = V 2 RT 


(1.8a) 


(1.8b) 


0 - 9 ) 


Associated with this relation (which is a simplified form of the equation 
obtained by taking the horizontal divergence of the equations of motion) is the 
"vorticity equation": 


= 

at 


k x • V (f+V 2 ty) + V • fVX + V • (fxk) ; 


t 




( 1 . 10 ) 


in which F represents the total dissipation forces, F r is the horizontal fric- 
tional stress force per unit mass, and f-g is the horizontal diffusion. 

The continuity equation (conservation of mass) is 


a_ 

3p 


'dpi _ 3_ 
dt J 3P 


dP 

dt 


= -V 


= v 2 x 


( 1 . 11 ) 


The upper boundary condition at Z = will be that dp/dt vanishes there. Let 

us define 


X = 



3X 

3P 


( 1 . 12 ) 


Equation (1.10) can then be rewritten as 
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+ V • (Fxk) 


0.13) 


v 2 M 

3t 


- k x Vip • V(f+V 2 ip) - V • fv 


fax) 

3P 


If we use Z = -InP as the vertical coordinate, the appropriate vertical 
advection velocity is 


w _ dZ 1_ dP_ 

dt P dt 

The continuity equation (1.11) in terms of W is: 

V-Pv + 3(PW)/3Z = 0 

From (1.11), (1.12) and (1.14) we get 3 [PW - V 2 x] / 3P = 0, or 

PW = V 2 x 


(1.14) 


0 . 15 ) 


(1.16) 


Boundary conditions on W are that W vanishes at and that it is given 
by orographic upslope motion at the bottom: 


Z = Z to P : W = 0 


Z = 0:W ::£ 2W = X 


V 4 


Vh 


(1.17) 


where h is the orography and v 4 is kxVip at the first interior level for Here 


H 0 = = 7 km (1.18) 

is a constant. 

1 -► 

Friction will be represented by a vertical Austausch, Fr = - 3x/3z = 

P 

-g3x/3p. Thus V*^rxk = |p- [7 • (j^-TXk)]. In the interior regions of the model 

(but not at the ground), in terms of a vertical momentum diffusion coefficient, 

A 

K , we set x = pK m 3(kxVi|>)/3z , giving 
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f 


■»-t=9?*k) - V-[2^K m 5|i] 


Replacing p by p/RT and replacing g/RT by 1/H 0 we get 


V-P^rxk] 

Po 


3V z ijj 
3 Z 


At the ground , we can set x equal to 0.003 p 0 |v|v, with j v"[ a suitable mean 
anemometer speed (5 m/sec ) and the anemometer vector wind v equal to a rota- 
ted (a = 22.5 degrees) fraction (0.5) of kxVip at the lowest interior level at 
which ip is defined 


Vnd = t°- 003 P^U°- 5 )^ cos a ^ “ sin a ^interior ~1 

l 0 - 19 ) 

V-C- ^xO grnd = -|^[.003|7|(0.5)c O s a]l7»* 1nterior J 

For H 0 = 7 km, |v| = 5m sec -1 and cos a = 0.925, the coefficient here has the 
val ue 10 ^ sec ^ . 

The conventional quasi-geostrophic Taylor-Ekman theory (Charney and Elias- 
sen, Tel lus , 1949, Vol . 1, No. 2, P.38) gives a corresponding term ("Ekman 
pumping") of 




sin 2 a}V 2 4i 


( 1 . 20 ) 


4 2 1 4 1 

For K m = 6.2x10 cm sec and f = 10 sec , the coefficient in this derivation is 
2. OxlC^sec - "' . To summarize the horizontal friction stress term, we can write 
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V-frxk = fp( PF ) 


3P 

Z>0- F = — 2- P 
L r FF K 3Z 
0 


2=0: F.= - k D 7** int 


( 1 . 21 ) 


where kg refers to the "surface drag coefficient" in (1.19) or (1.20). At 
Z = Z t , F will vanish (no stress). 


We want to note here that the pressure dependence of K m as used in our 
model is the same as for Kg, the coefficient of vertical diffusivity of ozone, 
to be introduced later. We, thus, defer discussion of K m and Kg values to the 
subsection entitled "Diffusivity for Ozone" but the results can be seen in 
Table 2.2. 

Horizontal diffusion in the model is primarily introduced at the upper 
(mesospheric) levels to account for possible affects of large scale horizontal 
transportation resulting from energy supplied by gravity wave propagating from 
below. Parameterzation is accomplished through introduction of a coefficient 
of horizontal momentum diffusion, A^, such that 

F d = A h V 2 v (1.21.1) 


where then 


v • 



A H k 


VxV 2 v 


A H [v 4 ^ 



(1 . 21 . 2 ) 


Here "a" denotes the earth's radius and all differential operators are assumed 
to be in the spherical surface. In our parameterization, A^ will be a function 
of pressure only and will serve as the diffusion coefficient for momentum, tem- 
perature and ozone. Again, discussion of the values to be associated with A^ 
will be deferred to the subsection on "Diffusivity of Ozone." 
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The next physical statement is the thermodynamic law d (entropy) /dt = 
rate of heating 4 - temperature. For our perfect gas system this would be 

C p dt K )] = x 1 K = q- = 7 ( 1 - 22 ) 

P 

where q is the rate of heating per unit mass and T the temperature. In terms 
of T, this becomes 


~ = -(kxV^-Vx) • VT - w|j - kWT + A h 7 2 T + 3- 


0.23) 


in which the term A^V 2 T represents the horizontal subgrid seal diffusion of 
temperature. We will, however, use a simplified form of this, obtained by 
ignoring Vifj*VT and by replacing T in W 3 T/ 3 Z and kWT by T, where T is the 
hori zontal average : 


T = T(p,t) + T"(A t) 
r-n/2 r 2 ir 


T = 


1 


4ira' 


cos 4 >d<|> 
-tt/2 


TdX; 1' = 0 


0 


(1.24) 


[ This definition of (~) and ( )' will be applicable to any variable. ] This 
greatly simplifies the computations, and is reasonably accurate because ~.‘j » 7;< 
and 3T'/3Z + kT" is generally small compared to 3T/3Z + kT. The result is 


§£ = - kxViji • VT - W(g + <T) + A H v 2 T + q/ C p 


(1.25) 
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However, this simplification has the result that we can no longer interpret 
(1.25) as forecasting T, the horizontally averaged T; this is because the hori- 
zontal average of (1.25) gives simply 


3T 

3t 


q/C 


P 


whereas the horizontal average of the exact equation (1.23) gives 


|T = 2_ . .w-T- - 


(1.26) 


showing the effect of vertical transports of entropy by the motion. We expect 
little change in T from the observed annual average T(Z), however, either with 
season or with changes in the ozone chemistry. [The effect of the latter will 
be discussed separately.] 

In passing, we note that 


3T 

3Z 


+ 


kT 


BI fll + sJ 
g c p J 


= % [£n(Tp _< )] 


Ni 

R 



(1.27) 


where N is the buoyancy frequency. 

Finally, we describe the basic form of the equation for the ( number density ) 
mixing ratio of a trace substance such as O^. Define 
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where n . is the number density of the i-th trace substance, n^ is the total num- 
ber density, assumed to be equivalent to the "normal" constituents and 

CO 2 since n^ is very small. 


"m S P/ kT 

-26 - 1 
k = Boltzman constant = 1.380x10 kilojoules deg" 


0 . 29 ) 


The equation for dy^/dt (the rate of change following the motion) is 


d Xi 3Xi - 9 Xi 

dt 1 = 5T * ( kxv *- v x) ' V Xl + V- 


1_ 

n 


fdn^l 


m 


dt 


+ 13 


p 3z 


pK 


d 3z 


+ 


where (dn.dt) is the net rate of local photochemical generation of the sub- 
stance (number per unit volume per unit time) and and A^ are the vertical 

p 

and horizontal eddy-diffusion coefficients [with dimensions (length) r time], 
respectively. and A^ will vary only with P. 

The vertical diffusion term can be rewritten by using the hydrostatic equa- 
tion as 



fs£ 

RT 


2 !^i] = 1_[_ 

3P J 3P l 



3Xi 

3Z 


•] 


(1.30) 


where we have again absorbed the variation of density with T into H 0 on the 
recognition that itself is not a precisely known quantity. (and the 
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momentum Austausch K m ) will be prescribed functions of P. The equation for x^ 
is now 


3x i 

-j^r- = -(kxV4)-Vx) * V X . 


3x i 

u 1 + 

W 3Z - 


1 

[ dn il 

n m 

[dt J 


|p [ " HT P 5T ] + A H v2 '-<i n -31 )- 


or 


9 Xi -1 ^ 3(PWx,) 

It = ’ f^ V, ( Pvx i) + 3F“ ] 


1 

Kl 

n m 

dt 


3x i 


0.32) 


3P 


C- HF p ST* + V 2 *l 


[having made use of (1.4) and (1.15) to obtain the last form]. 

The rate of change of x-j (the horizontal average) is obtained from the 
horizontal average of (1.32): 


3Xi_ 

at 


i 

dn . 

i 

l m 

dt 

\ ^ 


* U - 


K d p *1. 

Hf P TZ J 


0.33) 


The rate of change of x^ will, however, be obtained from a simplified form 
of (1.31), much as was done in the thermodynamic equation (1.25): 


3Xi - 3Xi 

__ = _ kxV^ • 7 X r - W^- 


1.34 


1 dn i 

+ [J Ll 

L n dt J c 
m 


K 


+ — f- — 1 p — 1 + A V 2 yT 
3P L Hi H 3Z J H V x l 


In contrast to T, where we are for the most part content to take T as given, we 
must predict X-j as we H as xj- Equation (1.33) will therefore be used as well 
as (1.34). 
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Presumably (1.33) need not be applied every time step in the numerical 
integration, y. being a slowly changing function of time. However, the term 
W xj must be put equal to zero at P = 1 to ensure no net creation of Xj by the 
large scale motion. 

The form of (dx^/dt) c is discussed later. However, a special treatment 
must be used for the lower boundary condition on the vertical eddy flux. For 
the case of ozone, Galbally ( Quart. J. Roy. Meteor. Soc. , 1971, P. 18) shows 
that in the very lowest layer the vertical flux (over land) is proportional to 
the ground concentration 


K IX « IX = d Y 
N d 3z H 0 3Z Xgrnd 


(1.35) 


(the surface destruction of ozone being proportional to 0^). The coefficient d 
has a value of about 1 cm sec""*. We will apply this formulation to the lowest 
layer in our model (0 < Z < AZ). Values of x are defined at the top of the 
layer (Z= AZ, j = J-l) and at the ground (Z = 0, j = J). Thus (dropping the 
i -subscript on x)-, 


f K j . 1 

J. lx 


K d] 

Ho 3Zj 

J-l/2 i 

H 0 AZ 


(x 

J-l/2 J-l 


-X ) = d X 
J J 


whence 


*J = Xj-1 


[1 + 


d H o AZ 

K d J J-l/2 


(1.36) 


0-37) 


and 


1-11 


(1.38) 


h 0 az 


J-l/2 


( Xj-TXj) = T 


d Xj-1 

+ (d HoAZ/K, 


Galbally cites values of the vertical number flux of ozone molecules over 

11 -2 1 

land in the range 1. to 6x10 mol cm sec . Aldaz ( J. Geo. Res. , 1969, P. 6943) 

11 -2 -1 

estimates a global average of 1 to 1.7x10 mol cm sec . Picking a represent- 

11 -2 -1 

ative value of 2x10 mol cm sec and equating this to n m K 3x/3z, we find, 
19-3 52 -1 

for n m = 4.55x10 cm and K = 10 cm see - , that a vertical gradient of ozone 
number mixing ratio of 


|| ~ 0.5xl0" 13 cm" 1 


5xl0~ 8 
10 km 


(1.39) 


is required. Galbally's data show a typical ground value for x of 5xl0 11 v4.5x 
19 -8 

10 ~ 10 . The typical inferred downward flux of ozone observed near the 

5 2 -1 

ground is compatible then with a tropospheric K of 10 cm sec and a tropopause 

(10 km) value for x of 6xl0~ 8 or a 10-km value for NO^ of (6xl0~ 8 ) x (8xl0 18 ) ~ 

10 -3 12 -3 

50x10 cm . This- value is not greatly inconsistent with values of 10 cm 

which seem characteristic of the tropopause level in the meridional cross sec- 
tion prepared by D. Wu from the data of Hering and Borden (1967). 

A special treatment of the minor species equation will be necessary at cer- 
tain levels. As an example, Lindzen and Goody ( J. Atmos. Sci . , 1965, P. 341) 
show that the photodissociation of ozone is extremely rapid at heights above ~ 45 
km, with a time constant becoming less than 1 hour. (They presumably use typi- 
cal values of incident solar radiation). The conventional methods of "time- 
stepping) equations such as (1.34) require a computational time step no longer 
than the characteristic physical times associated with terms on the right side 
of (1.34). Since the advective time scale is of the order of an hour or so, we 
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must consider replacing (1.33) and (1.34) at these levels by the equi 1 i bri um 
condition. 


dn 


Xi (X-j ) e q U i 1 < > dt ^ 


(1.40) 


For use in radiation computations, we need N.. , the number of molecules of 
X.j in the vertical column of unit cross section above a given pressure surface: 


N i “ 


n^dz = 


Xi" m dz ' Xj F P d2 


rP 

gk 


X-j d P 


where R = 287 kj ton""* deg ^ is the g.as constant for air. 
This gives numerically 


N. = 2.12x10 


29 


= 2.12x10 


25 


X^dP in (meter) 

5 _ 2 
X-jdP in (cm) 


-2 


(1.41) 


In the case of molecular oxygen, x-j I s taken as uniform and equal to 0.2096, 
gi vi ng 

N = 0. 4444x1 0^ 5 P cm' 2 (1.42) 

°2 
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The lower boundary condition for ozone 


The surface destruction process for ozone must give rise to a boundary 
layer effect in which the ozone changes rapidly from the free air value to a 
lower value at the ground. Our use of (approximately) 3-km height increments 
will not represent this adequately. Secondly, land and water surfaces differ 
markedly in their effect on surface ozone. Fortunately, it is possible to cor- 
rect for both of these complications by using the detailed analysis by P. Fabian 
and C. Junge (Global rate of ozone destruction at the earth's surface: Arch. 

Meteor. , Geo. , Biokl imat. , (a)-Meteor. u. Geo. 19, 161-72, 1970). The important 
point is to obtain the correct global ratio between the lower tropospheric value 
of ozone and the surface destruction rate, since the former may affect the ozone 
amount higher up and thereby the photochemical destruction rates. 

Fabian and Junge model the presumed boundary-layer ozone profile by stan- 
dard methods and make allowance for different properties of land, vegetation and 

water (and their global distribution) and for different wind speeds. They ar- 

1 0 2 

rive at a global surface destruction rate ranging from 3.1 to 5.6x10 mol /cm 

sec., the variation being due to uncertainty in choice of surface wind speed. 

1 1 3 

Using an average lower tropospheric value of 5x10 mol/cm for 0^, we have a 
global ratio of 


j _ 3.1 to 5.6x10 
d l = TT~ 


10 


- 0.08 cm sec 


-1 


5x10 


[Note that "measured" values of d at the ground range from .04 cm sec over 
water to land values of 0.6 (Aldaz) and 1.0 (Galbally) cm sec \] 

Referring to our model equations (1.37) and (1.38) 


1-14 



i H o AZ 


( xj-rxj) = — 
0 - 1/2 d 1+ 


%-! - H 

; - d x i 

fdHoAZ) 

Kr> 


we recognize the left side as the downward diffusive flux at the bottom of our 
model, (which must equal the surface destruction rate) and Xj _q as number 
density (mixing ratio) in the model corresponding most closely to the 5x10 
number density for the free air referred to above. Our model will not include 
different types of surface with their differing abilities to destroy ozone. 

This is alright since, because of the strong horizontal advection in the atmo- 
sphere, these differing surface properties affect primarily the local boundary 
layer profile and surface value rather than the local free troposphere values. 

We must use a correct global effect, however, and we get this by simply choosing 
a single model value for d such that the ratio of the destruction rate to the 
free air value in the model matches the global observed ratio, d. 


- = d, = 0.08 cm/sec 

, dH 0 AZ 1 

k d 

O I - ] 

For H 0 AZ = 3 Km, and K d = 10 m sec , this gives d = 0.105 cm sec . 
Fabian and Junge also discuss the ratio (their e) of the surface value of 0^ to 
the free air value, and obtain typical values of 0.35 over land and 0.85 over 
water. Our model now implies a single value for this ratio of 


JSl 

X J-1 



very comfortably located in the range inferred by Fabian and Junge. 
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Diffusivity for ozone 


Vertical diffusivity 

An upper estimate for K d in the troposphere and lower stratosphere can be 
obtained from measured ozone profiles by equating 


(48) 

29 




-12 2 

constant = 3.5x10 gm/cm /sec 


-1 2 

where y is the number density mixing ratio, and 3.5x10 is the product of the 

-2*3 

mass of an ozone molecule (7.9x10 gm) with the average global surface destruc- 
tion rate of 4.4x10^ mol cnf^sec"^ found by Fabian and Junge (reference cited 
earlier). Values of 3y/3z can be obtained from the middle latitude synthesis by 
Krueger and Minzer (A proposed mid-latitude ozone model for the U.S. Standard 
Atmosphere. X-651 -73—72 , Goddard Space Flight Center). (Similar values are 
obtained from using the 3-year average ozone profiles for Bedford and Green Bay 
that has been analyzed by D. Wu). The results of this calculation are shown in 
Table 1.1. 

At higher elevations ozone begins to no longer act as an inert tracer. 

Here we refer to computations by S. Wofsy and M. McElroy (On vertical mixing in 

the upper stratosphere and lower mesosphere. J . Geo. Res. , 78 , 2619-2624, 1973). 

-1 /2 

These authors combine (a_), the suggestion by Lindzen that K might vary as p 
because the velocities in gravity waves - a likely mixing process - tend to in- 
crease in this manner with height, and (b_) , measurements by Ehhalt of methane 

concentration. Using a chemical model, they find that Ehhalt' s measurements at 

3 2 -1 

50 km fit but with a distribution having small values of 2x10 cm sec at 

- 1/2 

16-20 km increasing to 2x10 at 80 km as (n^) - ; . We can model this simply by 
noting that n^ is proportional to p for constant T, and that p in turn is pro- 
portional to exp(-Z). 
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Table 1.1: Upper estimate for K d in the troposphere and lower stratosphere 


(km) 

f 48 ] 

[29 X j 


29 az v ' 

ill 

CL 

» cm 
\i sec 

2 

5x1 0 -8 




3 



-13 

0.5x10 

0.90xl0“ 3 

7. 80x1 0 4 

4 


6 




5 



0.5 

0.74 

9.40xl0 4 

6 


7 




7 



1.5 

0.59 

4. 00x1 0 4 

8 

10 




9 



6 

0.47 

1 .20xl0 4 

10 

22 




11 



15 

0.36 

0.65x1 0 4 

12 

52 




13 



15 

0.26 

0.90xl0 4 

14 


82 




15 



30 

0.19 

0.61x1 0 4 

16 

142 




17 



63 

0.15 

0.37x1 0 4 

18 

267 




19 



82 

O.llxlO -3 

0 . 39x1 0 4 

20 

431x1 0" 8 
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However, we believe that the three-dimensional model waves can sufficiently 
handle the vertical "diffusion" in the upper stratosphere and thus we have chosen 
the following distribution of K which essentially ignores vertical diffusion pro- 
cesses in the stratosphere. 

(a) Z < Z 0 = 0.6 : K = K 0 = IxlO 5 cm 2 /sec 

(b) Z 0 < Z < Z 10 = 2.4: . K r = IxlO 2 cm 2 /sec 

K = K° + [(Z-Zo)-(2-Z 10 -Z-Z 0 )] 1/2 

(c) Z > Z 10 : K = Ki 

The values used in the model as obtained from the formulas are shown in 
Table 1.2 along with values inferred from observed CH^ (Wofsy and McElroy) 
and ozone profiles. 
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Table 1.2: Values of K, from the model "formula" and inferred from 

a 2 

0 3 and CH^ profiles. K d units are cm /sec. 


z(km) 

Z 

Formula 

Inferred (0^) 

Inferred (CH^) 

3 

.357 

100. 0x1 0 3 

78.0xl0 3 

300. Oxl 0 3 

5 

.616 

90.0 

94.0 


7 

.892 

46.2 

40.0 


9 

1.171 

26.5 

12.0 


n 

1.492 

13.8 

6.5 


13 

1.802 

6.3 

9.0 


15 

2.120 

1.5 

6.1 


17 

2.433 

0.1 

3.7 

2.0xl0 3 

19 

2.719 

0.1 

3.9xl0 3 


20 

2.861 

0.1 


6.8 

30 

4,289 

0.1 


45.0 

40 

5.717 

0.1 


110.0 

50 

7.145 

0.1 


260.0 

60 

8.573 

0.1 


490.0 

70 

10.000 

O.lxlO 3 


940.0xl0 3 
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Horizontal diffusivity 


The purpose of the horizontal diffusion term in our model is principally 
to account for the energy conversion of vertically propagating gravity waves 
into horizontal transportation in the upper mesospheric levels of the model 
atmosphere and to act as a "sponge" to prohibit wave refelctions from the fixed 
(i.e., vertical motion = 0) upper boundary. Thus, we have chosen to specify A H 
as a function of 1 / 1 ^^ in the form 


A H - A const < Z / Z top> 6 C™ 25 ' 1 } 


(1.21.3: 


7 2-1 

in which A const is a specified constant. A value of A const = 6x10 [ms ] 
seems to work very well with the truncation specification of the current model 
version (18 wave) and has thus been adopted for all of its horizontal diffusion 
terms (i.e., for momentum, temperature and ozone diffusion). 

Profiles for both the vertical diffusion coefficients and for the hori- 
zontal diffusion coefficients as represented at the model levels are shown 
in Table 2.2. 
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2. Choice of vertical levels. 


We want the vertical domain to extend well above the actual ozone layer. 

We also want it high enough that there is some opportunity for the damping ef- 
fects of ozone and radiation to absorb mechanical energy generated in the baro- 
clinic processes of the lower atmosphere. On the other hand, we cannot for prac- 
tical reasons get involved in the more complicated processes of the lower -thermo- 
sphere. An upper limit of about 85 km seems reasonable. 

We obtain equal intervals in Z = -ZnP (P = pressure t 100 cb) by defining 

Z, =AZ(J-j) I 

J } i ■ 1, 2. 0. (2.1) 

p. = e -AZ(J-j)J 


j = 1 is at the "top" of our model atmosphere, and j = J at the bottom, whence 


Z 




A convenient choice is obtained by choosing 

e AZ = r, r = 3/2 

AZ = Inr = 0.40547 (2.2) 


so that 


Zi = 
Pi = 


Z top = (J-lttnr 


(2.3) 
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Successive pressure levels are separated by (roughly) 85/31 ~ 2.7 km. The 
relations 



are useful. At these levels, the following basic variables will be represented 

j = 1, 2, J : T . , W. , (x-)-; together with the heating rate, the photochemi- 

3 J ' J 

cal term, and the vertical turbulent fluxes of momentum. At the intermediate 
levels the streamf unction ip ■ will be represented 

J 



For convenience in notation, however, i p will be labeled with an interger sub- 
script according to the convention 


" P j+l/2^ ~ 


This results in the scheme as seen in Figure 2.1. 

Table 2.1 lists the values of the more basic variables for the choice 
r = 3/2, J = 32. Values of' 7 above 30 km were taken from the 1965 CIRA annual 
mean, values at lower elevation coming from data based on statistics gathered by 
the Planetary Circulation Project at M.I.T. (To be precise, they were obtained 
from the latter as shown in a figure based on them in the thesis by A. Hollings- 
worth). The static stability parameter S is defined later in equation (3.20). 

Table 2.2 contains the values for the vertical (K^, K^) and horizontal (A^) dif- 

o 

fusion coefficients (in units of m /s) at the same model levels as shown in 
Table 2.1. The basis for these values is discussed in the previous section. 
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Figure 2.1 


////////////// 


////////////// 


Vertical levels of the model and the location on these 
levels of the model variables. 


j 


1 

P 1 

z i 

W 1 c=0) 


tx,), 

T 1 

F 1 





*1 

^1 




G 1 

2 

P 2 

Z 2 

W 2 


(x f ) 2 

T 2 

F 2 





i>2 

^2 




G 2 

3 

P 3 

Z 3 

W 3 


(Xi), 

T 3 

F 3 



j-1 

P 0-l 

Z J-1 

W 3-l 

^ X i ^j-1 

T 3-l 

F j-1 




*3-1 



G j-1 

j 

P 3 

Z 3 

W 3 

(Xi)j 

T J 

F 3 




*j 



G 3 

j+1 

P j + 1 

Z j + 1 

Vl 

• 

(Xi ) J+1 

T j + 1 

F j + 1 




• 





J_1 P J-1 Z J-1 W J-1 


(x i ^J-l T J-1 F J-1 

^J-l 

^J-1 

G J- 

J 'j Z J W J 


(Xi)j T J F J 
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TABLE 2.1: Vertical layer designations and mean vertical profiles 


j 

+ 

p j 

z j 



n (cm) 3 
m v ' 

s j 

1 

.00000348 

12.569 

85.4 

181 .0 

1 39x1 0 1 2 

173 

2 

.00000522 

12.164 

83.2 

181.0 

209 

173 

3 

.00000782 

11.758 

81.0 

181.0 

313 

162 

4 

.0000117 

11.353 

78.8 

183.5 

463 

130 

5 

.0000176 

10.948 

76.5 

192.0 

664 

111 

6 

.0000264 

10.542 

74.1 

201 .0 

952x1 0 1 2 

114 

7 

.0000396 

10.137 

71.6 

211.0 

1 36x1 0^ 3 

137x10 

8 

.0000594 

9.731 

69.0 

219.0 

196 

144 

9 

.0000891 

9.326 

66.3 

226.5 

285 

154 

10 

.000134 

8.920 

63.5 

234.0 

415 

161 

11 

.000200 

8.515 

60.6 

241.5 

600 

166 

12 

.000301 

8.109 

57.6 

249.6 

877x1 0 1 3 

167 

13 

.000451 

7.704 

54.5 

258.5 

126xl0 14 

174 

14 

.000677 

7.298 

51.4 

267.0 

184 

217 

15 

.00101 

6.893 

48.2 

267.5 

274 

277 

16 

.00152 

6.488 

45.0 

261.5 

421 

302 

17 

.00228 

6.082 

41.9 

254.5 

649x1 0 1 4 

295 

18 

.00343 

5.677 

38.8 

248.5 

lOOxlO 1 5 

285 

19 

.00514 

5.271 

35.9 

242.5 

154 

277 

20 

.00771 

4.866 

33.0 

237.0 

236 

272 

21 

.0116 

4.460 

30.2 

231 .0 

364 

269 

22 

.0173 

4.055 

27.5 

225.0 

557 

261 

23 

.0260 

3.649 

24.8 

219.5 

855x1 0 1 5 

251 

24 

.0390 

3.244 

22.2 

214.5 

132xl0 lb 

237 

25 

.0585 

2.838 

19.6 

211.5 

201 

217 

26 

.0878 

2.433 

17.1 

210.5 

302 

iga 

27 

.132 

2.027 

14.6 

213.0 

499 

155 

28 

.198 

1.622 

12.0 

222.0 

646 

125 

29 

.296 

1.216 

9.3 

234.0 

91 3x1 0 4 6 

116 

30 

.444 

0.811 

6.4 

253.0 

1 30x1 O 1 7 

104 

31 

.667 

0.405 

3.4 

272.0 

182 

105 

32 

1.000 

0.0 

0.1 

287.0 

265x1 O 1 7 

122 
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TABLE 2.2: Profiles for the vertical and horizontal diffusion coefficients 


j 

z j 

Vo 

(m 2 /s) 

a h 

(m 2 /s) 

1 

12.569 

9.997x10° 

6. 000x1 0 7 

2 

12.164 

8.044 

4.928 

3 

11.758 

6.449 

4.021 

4 

11.353 

5.147 

3.258 

5 

10.948 

4.033 

2.619 

6 

10.542 

3.201 

2.088 

7 

10.137 

2.506 

1.650 

8 

9.731 

1.928 

1.292 

9 

9.326 

1.455 

1.001 

10 

8.920 

1.069 

7.665xl0 6 

11 

8.515 

7. 544x1 0" 1 

5.798 

12 

8.109 

4.972 

4.327 

13 

7.704 

2.872 

3.108 

14 

7.298 

1.157 

2.299 

15 

6.893 

1. 000x1 0 -2 

1.632 

16 

6.488 

1.000 

1.134 

17 

6.082 

1.000 

7.701x1 0 5 

18 

5.677 

1.000 

5.090 

19 

5.271 

1.000 

3.263 

20 

4.866 

1.000 

2.019 

21 

4.460 

1.000 

1.198 

22 

4.055 

1.000 

6.760xl0 4 

23 

3.649 

1.000 

3.593 

24 

3.244 

1.000 

1.772 

25 

2.838 

1.000 

7.954x10° 

26 

2.433 

1 .000x1 o' 2 

3.154 

27 

2.027 

2.265x1 0 _1 

1.056 

28 

1.622 

9.91 7x1 0" 1 

2. 769xl0 2 

29 

1.216 

2.473 

4.928X10 1 

30 

0.811 

5.308 

4.327 

31 

0.405 

10.000 

6.760xl0- 2 
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3. Non-dimensional finite-difference equations 


In this section we write the basic equations in a non-dimensional form (pri- 
marily to simplify the dynamical computations) and simultaneously introduce the 
vertical finite-difference representation defined in Section 2. We define 

y = sin<j) 


V ( d i m ) 

V 2 (dim) 
i|;(dim) 


— V(non-dim) 

d 


V 2 (non-dim) 


a^ 

29a z ^(non-dim) 
2£2a 2 X(non-dim) 


X(dim) 
t(dim) = ^ t(non-dim) 


1 


(3.1) 


W(dim) = 29 W(non-dim) 

Aj_j ( dim) = 2Qa 2 Aj_|(non-dim) 


T(dim) = (4fi 2 a 2 /R) T(non-dim) + (4ft 2 a 2 /R) T(non-dim) 


In the last expression T (dim) is the "total" temperature in absolute degrees, 

T = T(Z) is the "standard atmosphere" temperature (also in degrees) given in the 
table at the end of Section 2, while the quantity (4ft 2 a 2 /R) T (non-dim) is the 
(deviation from the horizontal mean) variable T appearing in (1.25), having a 
zero horizontal average. [The total T (dim) is, of course, used in all chemical 
computations] . 

9 = 2tt/8. 64x10^ rad sec 
a = 6.371x1 0^ meters 
R - 287 kj ton' 1 deg' 1 
Cp = ( 7/2 ) R 
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One day, (2-rr/Q) secs, corresponds to 


At(non-dim) = = 4ir 


(3.3) 


The non-dimensional V 2 operator is 


o2/ \ - } _ 3 2 ( ) + _J i_ f cos(b ^ - H 

V ' ' C0S 2 <f> dX 7 ^ COS(f> 3<j> Lcos * 34> -* 


(3.4) 


The relation 


PW = V 2 X (1.16) 

between W and X can be used to eliminate X in favor of W [in equation (1.13)] 
by defining the inverse Laplacian operator 


L = V“ 2 
X = PLW 


(3.5) 


We also have 

C = V 2 lj) ; 4> = Lz 


(3.6) 


A further convenient arrangement is useful for evaluating terms of the form 
3(PF)/3P, which appears in the vertical diffusion terms for vorticity and trace 
substances and in the term 

fjr = fpCP(LW)] 


in the vorticity equation (1.13). We have 
P . 


[^(PF)] j 


j+1/2 


F j + 1 /2 " P ,j-l/2 F j-l/2 


P j+l/2 " P j - 1 / 2 


^r-l )F j+l/2 ' ^r-1 


)F j-l/2 


(3.7) 


where we have made use of (2.4). 
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The horizontal advection of a quantity F can be written as the Jacobian 


-v^ • VF = - kx7^ • VF = 


3F 3^ 
3A 3y 


djjJ 

3 a 


3F 

3y 




(3.8) 


The non-dimensional form of the vorticity equation (1.13), with regard to 
the subscript labelling defined in Section 2, together with equation (1.21) and 
(3.5) - (3.8) is as follows: 


For j = 1 , 2 , . . . , 0-1 : 


Jfu+Cj.tj) - V-ipVLt^Wj^ - (^Tj-lWj] + 

- (r?r)F J+1 - trJrifj + 

F, = 0 
F J = 

F j ■ (J = 2. 3. .... J-l) 

Ej * (K m )j t [HS20AZ] 

D = k Q t 2Q 

w 1 = 0 

W J = ^J-l^ 

0 


(3.9) 

(3.10) 

(3.11) 

(3.12) 

(3.13) 
(3-14) 

(3.15) 

(3.16) 
(3.17) 


The non-dimensional form of the "thermal wind equation" (1.9) becomes for 
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j = 2, 3, ..., J-l: 

= -V 2 TjAZ ( 3 . 18 ) 

The non-dimensional form of the thermal equation (1.25) becomes for 


j = 2, 3, . . . , J-l : 


Hi 

3t 


J ‘V - S J W 0 + V/‘ t j + [ o§ 


T^rl 


where 




(3.19) 


(3.20) 


is tabulated at the end of Section 2. Note that q., the rate of heating per 

J 

unit mass, is still in dimensional form in (3.19). It is considered later in 
Section F. 

The trace substance is, for 


j = 0 o , j Q +l , , J-l : 


3X 


3t 


1 - j j(xj. Kj-tj.,) - Vaz> + ( ^t )g j ' ( r^r )G j-i + Vj^xj 


m 


(3.21) 


WXj+lV ; for j-j 0 . .... J-2 


Dj - (K d ) j+ , /2 4- (2 «h;aZ) 


(3.22) 


[The vertical diffusion coefficient K. is defined at the Z.-levels corresponding 

d J 
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to j = integer plus 1/2, whereas the vertical exchange coefficient K m for vor- 
ticity, appearing in (3.14), is defined at interger values of j.] At the bottom, 
the relation (1.38) gives 


G 


J-l 


Xj-1 


2QHo + 


2QHqAZ ' 
(K d } J-l/2 


The integer j sets the level above which (3.21) may be replaced by 
ical equilibrium statement, as discussed near the end of Section 1. 


(3.23) 


a photochem- 
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4. Model chemistry 


i cal 
1 . 
2 . 

3. 

4. 

5. 

6 . 

7. 

8 . 

9. 

10 . 


The distribution of ozone in the model is determined by the following chem- 
reactions (with rate constants given by Baulch et al_. , 1982) 

1560 

NO + 0 3 -* N0 2 + 0 2 k 3 = 3.6 x 10" 12 e T 


N0 2 + 0 - NO + 0 2 


k 4 = 9.3 x 10' 


•12 


N0 2 + hv ■* NO + 0 


J 


no 2 


OH + 0 3 -> H0 2 + 0 2 
OH + 0 -*■ H + 0 2 
H + 0 2 + M -> H0 2 + M 
H + 0 3 -*■ OH + 0 2 
ho 2 + 0 -*■ OH + o 2 
H0 2 + 0 3 - OH + 20 2 
0 3 + 0 -*■ 20 2 


1000 

-12 j 

k ]2 = 1.9 x 10 u e 

no 

k ]4 = 2.3 x 10" 11 e T 


k ]6 = 5.9 x 10 


-32 


f -r- ^ 


300 


'17 


480 

1 .4 x 10" 10 e ' T 


k^g = 3.7 x 10 


■11 


600 

-14 - T 

k 2Q = 1.4 x 10 e 


2300 

k 22 = 1.8 x 10" 11 e ‘ T 
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11. + hv -*■ 0 + 0 2 

12. 0 + 0 2 + M->0 3 + M 

13. C >2 + hv -*■ 20 

14. CIO + 0 - Cl + 0 2 

15. CIO + NO - Cl + N0 2 

18. 0 + 0 + M -»■ 0 2 + M 

19. NO + H0 2 + N0 2 + OH 


C 25 = 6.2 x 10 

S 

k 31 = 7.5 x 10' 11 
k 32 = 6.2 x 10“ 12 
C 37 = 4.8 x 10' 38 
k 33 = 3.7 x 10 12 



- 2.0 


e 


120 

T 


e 


294 

T 


e 


900 

T 


e 


240 

T 


We assume chemical equilibrium determines the N0/N0 2 balance, the Cl /CIO 


balance, and the 0H/0H 2 /H balance. Thus, we assume 

k 3 [N0][0 3 ] + k 32 [C10][NO] + k 38 [NO][H0 2 ] = k 4 [N0 2 J[0] + [NO^J^ (4.1) 
k 14 [0H][0] = k 16 [H][0 2 ][M] + k 1? CH][0 3 ] (4.2) 

k 17 [H][0 3 ] + k 18 [H0 2 ][0] + k 38 [N0][HO 2 ] + k 2Q [H0 2 ][0 3 ] = k 12 [0H][0 3 ] + 

k 14 [0H][0]. (4.3) 

We also assume that the 0/0 3 balance is determined by the principal terms only: 

[O3] Jq 3 = -€25 Co] [O2] Cm 3 ( 4 - 4 ) 
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With these assumptions the chemical loss of odd oxygen is given by 


CM] 



= 2[0 2 ]Jq 2 - 2k 4 [N0 2 ][0] - 2k 12 [0H][0 3 ] - 2 fc [0H][0 3 ] 

- 2k u [0H][0] - 2k 22 [0 3 ][0] - 2k 31 [C10][0] 

- 2£ 37 [0] 2 [M] + 2k 38 [N0][H0 2 ]. 


(4.5) 


Two-dimensional daily-averaged solsticial distributions of the mixing ratios 
of N0 X , CIO, and OH have been provided by Dak Sze (see Figures 1-3). We have as- 
sumed that these distributions are appropriate for January 1 and that the distri- 
butions evolve sinusoidally in an annual wave throughout the year. Most of the 
chemistry in the model is active only during the hours of sunlight. Therefore, 
we solve the above chemical equilibrium equations for the average species concen- 
trations during daylight hours. The equations result in a quadratic equation for 
[NO]. The solution to this equation then readily provides daytime values of 
[N0 2 ] and [H0 2 ]. These daytime values are then inserted into equation (4.5) and 
the entire right hand side of the equation is multiplied by h s /ir (see Section 5). 

There exist conditions for which the above determination of the chemical 
loss of odd oxygen may be inappropriate. Note, in particular, that we are assum- 
ing the complete absence of chemistry in the polar night. We have ignored the 
photodissociation of NO which occurs above 50 km; this means that the calculated 
N0/N0 2 ratio will be incorrect at high altitudes, but this should not affect the 
other ratios we calculate. 

The numerical model uses half-hour time steps (At). Stability of the model 
then requires a way to predict 0^ when its chemical time constant is less than 
or equal to approximately 0.5 hours. The following procedure is used: 


4-3 



At the start of each N-cycle we calculate 0.42[J n ] where [J n ] is the 

u 2 u 2 

zonal -mean daily-average photodissociation rate of molecular oxygen at each of 
the Gaussian latitudes. We also calculate 


TC = [ Xn ] - 0.42[J n ] • (0.7 x tt x At). (4.6) 

u x u 2 

At those latitudes and levels where TC > 0, a chemical equilibrium value of 0 

A 

is calculated by setting the RHS of equation (4.5) equal to zero and with the 
O/O-j ratio being determined by equation (4.4). At these locations (3 xq /3t ) c 

X 

is simultaneously set equal to zero. A prediction is then made on Xn in the 

u x 

spectral domain. Next, Xq is transformed back to the grid domain and at the 

u x 

locations where chemical equilibrium applies, the predicted value of x 0 is then 

u x 

partitioned into ozone and atomic oxygen according to equation (4.4). Finally, 

X Q and Xq are transformed back into the spectral domain. 
u x u 3 

The effect of the above procedure on the odd oxygen budget has been analyzed 
in Cunnold (1983) for the 6 wave model. In that report it was noted that the 
boundary between the regions where photochemical equilibrium is assumed and 
where it is not, constitutes a discontinuity which, because of the spectral 
representation of model variables, produces minor numerical effects on the ozone 
budget at all latitudes. We found, furthermore, that these effects are mini- 
mized if the chemical computation starts from the Xq 9 rid field and not from 

u x 

X n field. This is because, at high altitudes, strong gradients of Xn (but not 
' u 3 3 

of Xq ) exist across the polar night boundary which are poorly represented by 
x 

the spectral model . 

In treating the chemistry we use p/kT for the concentration of air molecules 
in the atmosphere and 0.21 times this value for the concentration of oxygen mole- 


k 
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cules . To increase the efficiency of the computation of reaction rates 
are calculated at every time step, we approximate each rate constant as 
duct of the variables exp(-125/T) and exp(-300/T). 


which 
a pro- 
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5. The calculation of photodissociation rates and the associated heatin 


Photodissociation rates and heating rates due to the absorption of solar 
radiation by molecular oxygen, ozone, and nitrogen dioxide are evaluated in the 
model . 

The vertically overhead column of ozone is first evaluated through the 
integration 


N n (P.) = 2.1156 x 10 25 J \ dP cm -2 . (5.1) 

u 3 J 0 0 3 


We assume that Xn varies linearly with Z in the interval j to j+1 . Then 
u 3 





J-l 

E 

0=1 




(5.2) 


To increase computational efficiency the longitudinal variation of the 
photodissociation rates is approximated using the following procedure. At each 
time step the mean and standard deviation of the column ozone is obtained at 
each latitude and level. The photodissociation and heating rates are then eval- 
uated for the two ozone columns equal to the mean ±1 standard deviation. For 
example, the photodissociation rate of ozone is 



a Q (A i )F(X i ) ( 


‘° 3 Ui)X °3 


(5.3) 


where the absorption coefficients of molecular oxygen (c^) and ozone (ciq^) and 
the integrated photon flux in each wavelength interval (F(X.)) are given in 
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Table 5.1. The tabulated values of F have been taken from page B-10 of "The 
Stratosphere, 1981." The ozone absorption coefficients are from Inn and Tanaka 
but adjusted by the Vigroux temperature correction factors for -44°C as recom- 
- mended by Klenk (1981). At visible wavelengths, the values are from Ackermann 
(1971). The absorption coefficients of NO ^ are from Table A-5 of "The Strato- 
sphere, 1981." The absorption cross sections of in the Hartley continuum is 
from Table 1-33 of "The Stratosphere: Present and Future, 1979." In the Schumann- 
Runge bands we use the mean of the high and low values given by Kockarts (1971), 
a choice which is consistent with the recent work of Frederick and Hudson (1980). 
To account for the backscattering of energy from the lower atmosphere, we in- 

o 

crease all of the solar fluxes at wavelengths greater than 31 25 A by a factor of 
1.5. 

Heating rates are computed in a similar manner. The heating rate, due to 
the absorption of ultraviolet and visible radiation by ozone is given by 


, + a o /V x o 9 

h- 9-453 * 10 7 x 0 „?<* 0,<V F <*1>X7 e ' L 


(5.4) 


This is in dimensionless units and is a term in equation (T.25). Here Xn ‘‘ s the 

u 3 

ozone volume mixing ratio. 

In equations (5.3) and (5.4), X n and X n are the column densities of ozone 

u 3 U 2 

and molecular oxygen. (For typical globally-averaged values, see Table 5.2). 


X 0 = ^0 sec ^ cm 

Xg = 2.1156 x 10 25 Pj sec ij; cm 
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TABLE 5.1: Flux of solar photons, q, at one All, absorption cross 

section of 0 2 and of 0 3 , a(0 2 ) and c(0 3 ), for wavelength intervals 
AA and wavenumber intervals Av. 


No. AA(A) Av(cm~. 1 ) q(cm 2 s" x ) a(02)(cm 2 ) a(0 3 )(cm 2 ) o(N0 2 )(cm 2 ) , 


45 

1 . 

754-1 . 

739 

57. 

000-57. 

500 

1 .33X10 11 

2 . 

74x10' 

1 9 





46 

1 . 

770-1. 

754 

56. 

500-57. 

000 

1.61 

1 . 

40x10" 

20 





47 

1 . 

786-1 . 

770 

56. 

000-56. 

500 

2.00 

9. 

00x10 " 

2 1 





48 

1 . 

802-1. 

736 

55. 

500-56. 

000 

2.20 

5. 

00 






49 

1 . 

818-1 . 

802 

55. 

000-55. 

500 

3.18 

3. 

20 






50 

1 . 

835-1. 

818 

54. 

500-55. 

000 

3.30 

2. 

00 






51 

1 . 

852-1 . 

835 

54. 

000-54. 

500 

3.14 

1 . 

20 






52 

1 . 

869-1 . 

852 

53. 

500-54. 

000 

3.71 

5. 

00x1 O' 

22 


19 



53 

1 . 

887-1. 

869 

53. 

000-53. 

500 

4.96 

3. 

00 


5.70x10" 



54 

1 . 

905-1. 

887 

52. 

500-53. 

000 

5.57 

1 . 

30 


5.15 




55 

1 . 

923-1 . 

905 

52. 

000-52. 

500 

6.34 

7. 

00x10 ' 

‘23 

4.66 




56 

1 . 

942-1. 

923 

51. 

500-52. 

000 

6.53 

4. 

50 


4.25 




57 

1 . 

961-1 . 

942 

51. 

000-51. 

500 

9.01 

2 . 

90 


3.99 




58 

1 . 

980-1 . 

961 

50. 

500-51. 

000 

1 . 02xl0 12 

2 . 

,10 


3.58 




59 

2 . 

000-1 . 

980 

50. 

000-50. 

500 

1 .15 

1 . 

,70 


3.20 




60 

2 . 

020 - 2 . 

000 

49. 

500-50. 

000 

1.40 

1 

.50 


3.09 



25xl0' 19 

61 

2 . 

041-2. 

020 

49. 

000-49. 

500 

1 .69 

1 

.25 


3.09 


3. 

62 

2 . 

062-2. 

,041 

48. 

500-49. 

000 

2.07 

1 

.00 


3.35 


3. 

75 

63 

2 . 

083-2. 

.062 

43. 

,000-48. 

,500 

2.52 

9 

.80x10" 

'24 

4.04 


3. 

79 

64 

2 . 

.105-2. 

,083 

47. 

,500-48. 

.000 

4.21 

9 

.20 


4.91 


3. 

85 

65 

2 . 

,128-2. 

,105 

47, 

,000-47. 

,500 

7.23 

8 

.50 


6.37 


3. 

92 

66 

2 . 

.150-2, 

.128 

46. 

.500-47. 

.000 

7.79 

7 

.85 


8.32 

’ 1 8 

3. 

98 

67 

2 . 

.174-2. 

.150 

46. 

.000-46, 

.500 

8.45 

7 

.05 


1 .09x10' 

4. 

,01 

68 

2 . 

.198-2. 

.174 

45, 

.500-46, 

.000 

1 . 05x1 0 1 3 

6 

.15 


1.46 


3. 

,98 

69 

2 , 

. 222 - 2 , 

.193 

45, 

.000-45, 

.500 

1.19 

5 

.50 


1 .88 


3. 

,81 

70 

2 , 

.247-2 

.222 

44, 

.500-45, 

.000 

1.51 

4 

.75 


2.41 


3. 

.46 

71 

2 

.273-2 

.247 

44 

.000-44 

.500 

1 .33 

4 

.05 


3.08 


3, 

.08 

72 

2 

.299-2 

.273 

43 

.500-44 

.000 

1.31 

3 

.35 


3.86 


2 , 

.67 

73 

2 

.326-2 

.299 

43, 

.000-43 

.500 

1.51 

2 

.70 


4.71 


2 . 

.33 

74 

2 

.353-2 

.326 

42 

.500-43 

.000 

1.32 

2 

.20 


5.62 


1 

.68 

75 

2 

.381-2 

.353 

42 

.000-42 

.500 

1.50 

1 

.65 


6.72 


1 

.21 

76 

2 

.410-2 

.381 

41 

.500-42 

.000 

1.34 

1 

.20 


7.62 


7 

.50xl0 _2C 

77 

2 

.439-2 

.410 

41 

.000-41 

.500 

2.02 

7 

.50x10 

“2 5 

8.66 


0 

.22 

78 

2 

.469-2 

.439 

40 

.500-41 

.000 

1 .82 




9.74 

" 1 7 

4 

.20 

79 

2 

.500-2 

.469 

40 

.000-40 

.500 

1.88 




1 .04x10' 

3 

.30 

80 

2 

.532-2 

.500 

39 

.500-40 

.000 

1 .83 




1.09 


2 

.36 

81 

2 

.564-2 

.532 

39 

.000-39 

.500 

2.25 




1.10 


1 

.45 

82 

2 

.597-2 

.564 

38 

.500-39 

.000 

4.65 




1 .07 


1 

.72 

83 

2 

.632-2 

.597 

38 

.000-33 

.500 

4.44 




1.02 


1 

.95 

84 

2 

.667-2 

.632 

37 

.500-38 

.000 

1 ,07xl0 14 




9.40x10 

“ 1 3 

2 

.05 

85 

2 

.703-2 

.667 

37 

.000-37 

.500 

1.18 




8.24 


2 

.80 
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Table 5.1 Continued 


No. 

AA(A) 

Avtcm" 1 ) 

q(cm 2 s a(0 2 )(cm 2 ) 

c(0 3 ) (cm 2 ) 

a(N0 2 ) (cm 2 ) 

36 

2.740-2.703 

36.500-37.000 

l.OSxlO 14 

6.69xl0" 18 

3 . 50x1 0 _2 0 

87 

2.778-2.740 

36.000-36.500 

1.04 

5.13 

4.32 

88 

2.817-2.778 

35.500-36.000 

7 . 54x1 0 1 3 

3.75 

5.60 

89 

2.857-2.817 

35.000-35.500 

1 . 48x1 0 1 ^ 

3.09 

6.30 

90 

2.899-2.857 

34.500-35.000 

2.17 

1.62 

6.77 

91 

2.941-2.899 

34.000-34.500 

3.46 

1 .00 

7.50 

92 

2.985-2.941 

33.500-34.000 

3.39 

5.31xl0“ 19 

9.10 

93 

3.030-2.985 

33.000-33.500 

3.24 

3.85 

1.17x10 19 

94 

3.077-3.030 

32.500-33.000 

4.40 

1 .66 

1 .70 

95 

3. 100(±25) 

32.520-32.000 

4.95 

9 . 37x1 0~ 2 0 

1 .83 

96 

3.150 

32.000-31.496 

5.83 

4.52 

2.19 

97 

3.200 

31.496-31 .008 

6.22 

2.53 

2.35 

98 

3.250 

31.008-30.534 

6.96 

1.19 

2.54 

99 

3.300 

30.534-30.075 

8.61 

5.54xl0" 21 

2.91 

100 

3.350 

30.075-29.630 

8.15 

2.73 

3.14 

101 

3.400 

29.630-29.197 

8.94 

1.38 

3.23 

102 

3.450 

29.197-28.777 

8.44 

8. 39xl0” 22 

3.43 

103 

3. 725±250 

28.777-25.127 

9.85 

4.30xl0" 23 

4.95 

104 

4.225 

25.127-22.346 

1 .88x1 0 a 6 

5.95 

*5.30 

105 

4.725 

22.346-20.100 

2.40 

5.29x1 0" 22 

*3.30 

106 

5.225 

20.100-18.265 

2.46 

2.11x1 0" 2 1 

*2.00 

107 

5.725 

18.265-16.736 

2.61 

4.21 

*7.00x10" 2 ° 

108 

6.225 

16.736-15.444 

2.63 

3.80 


109 

6.900±425 

15.444-13.605 

4.32 

1 .28 



* Photodissociation of NC^ is assumed to occur for X < 3975A. 
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TABLE 5.2: Typical concentrations of ozone, molecular oxygen, and 

nitrogen dioxide. 


Level 

Height 

km 

n(0 2 ) 
(cm -3 ) 

Column Cone. 
No^om' 2 ) 

n(0 3 ) 

cm” 3 

Column Cone. 
N* 3 (cm' 2 ) 

n(N0 2 ) 

cm” 3 

1 

71 .6 

2.86xl0 14 

1 .76xl0 20 

2.8x10 s 

8. 12x1 0 1 3 


2 

69.0 

4.11 

2.64 

6.4 

1 .97x10 14 


3 

66.3 

5.99 

3.96 

1 .5xl0 9 

4.70 


4 

63.5 

8.72 

5.94 

3.5 

1 . 13x1 0 1 5 


5 

60.6 

1.26xl0 15 

8.91 

6.5 

2.52 


6 

57.6 

1.84 

1 . 34x1 0 2 1 

1.2xl0 10 

5.16 


7 

54.5 

2.65 

2.01 

2.3 

1 .03xl0 16 


8 

51 .4 

3.86 

3.01 

4.6 

2.08 


9 

43.2 

5.75 

4.51 

9.9 

4.30 


10 

45.0 

8.84 

6.77 

2 . 1 5x1 0 1 1 

9.01 

8. 84x1 0 6 

11 

41 .9 

1 . 36x1 0 1 6 

1.02xl0 22 

4.1 

1.84xl0 17 

7 . 79x1 0 7 

12 

38.8 

2.10 

1.52 

7.6 

3.55 

4.00x10 s 

13 

35.9 

3.23 

2.28 

1 . 3x1 0 1 2 

6.50 

1 .39x1 0 9 

14 

33.0 

4.96 

3.43 

1.8 

1 . 09x1 0 1 8 

2.60 

15 

30.2 

7.64 

5.14 

2.45 

1.68 

3.27 

16 

27.5 

1 . 1 7x1 0 1 7 

7.71 

3.4 

2.46 

3.40 

17 

24.8 

1.80 

1 .16xl0 23 

4.3 

3.48 

3.68 

18 

22.2 

2.84 

1 .73 

4.8 

4.66 

3.96 

19 

19.6 

4.22 

2.60 

4.65 

5.87 

4.40 

20 

17.1 

6.34 

3.90 

3.6 

6.92 

6.04 

21 

14.6 

9.43 

5.85 

2.5 

7.68 

8.08 

22 

12.0 

1 . 34x1 0 1 a 

8.78 

2.0 

8.25 

1 . 1 5x10 1 0 

23 

9.3 

1.92 

1 . 32xl0 2i+ 

1 .0 

8.66 

1.72 

24 

6.4 

2.73 

1 .98 

e.oxio 11 

8.88 

2.58 

25 

3.4 

3.82 

2.96 

6.0 

9.05 

3.87 

26 

0.1 

5.56 

4.44 

6.0 

9.25 

5.80 
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if; is the solar zenith angle which is determined from 


j cos if; = sin<f> sind + cos<f> cos6 cos h (5.5) 

where <f> is latitude, d is the solar declination and h is the hour angle (mea- 
sured from local noon). 

An additional approximation is used in deriving photodissociation and 
heating rates. The equations require average photodissociation rates during 
daylight hours and heating rates averaged over a 24-hour period. This averaging 
is accomplished by evaluating the rates at two preselected hour angles. If h $ 
is the hour angle of sunset, it is given by 


cos h = - tan<j> tan<5. 

Considering only the Northern Hemisphere (tan <p > 0), we have 


tan<{> tand < - 1 : 
-1 < tan<j> tand < 0 : 

0 < tanc); tand < 1 : 

1 < tan<j> tand 


h $ e 0 (polar night) 
0 < h $ < ^ 

\ < h s < TT 


h $ = tt (polar day) 


The daylight-average photodissociation and heating rates are determined as the 
average of the values at h = h^/4 and 3h s /4. Empirically, we have found that by 
dividing the result by 1.025 we remove a slight bias in the results. Our studies 
indicate, moreover, that the approximation is accurate to better than 5% (see 
also, Cogley and Borucki , 1976). Twenty-four hour averages are determined by 
multiplying by h^/iT . 
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Using the above procedure, the photodissociation and heating rates are de- 
termined at each latitude and level for ozone columns equal to the zonal mean 
([X]) ± 1 longitudinal standard deviation (a). Suppose the results are J + and 
J_. The photodissociation and heating rates produced by ozone and molecular 
oxygen are calculated at each longitude from 


J = 




X-[X] 

2a 


V, J 


( 5 . 6 ) 


NC >2 heating and dissociation rates are considered to depend little on the 0^ 

1 /? 

column. For these terms (J + J_) is used at all longitudes. 
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6 . Tropospheric and infra-red stratospheric heating rates 


In the troposphere and lower stratosphere, we follow Trenberth (Ph.D. 
thesis, M.I.T., 1972) in setting, at each level j, 


II 

at 


+ V T j 


- T) 


( 6 . 1 ) 


* 


where h- is a "Newtonian" cooling coefficient, T. is a hypothetical 
J J 

temperature field, and T is the temperature predicted by the model. 

the above equation have a zero horizontal average.) For the zonally 

* 


part, Trenberth divided T into an annual average term (symmetric in 
and a seasonal term (an odd function of latitude): 


equi 1 ibri urn 
(All T's in 
symmetri c 
latitude) 


T j*(u.t) = Tj2 P°(u) + C T j* p i(p> + T j3 ( , 3(p)]sin[|f 5 - (t-f,)]. 

Here y is the sine of the latitude and P^ are the Legendre functions 
so that 


( 6 . 2 ) 
normal i zed 


f 1 2 

P^dy = 2 

J -1 

P°(u) = 

p°(u) - ”1 ( v ! - 1) 

P°(y) = ^ (5y 2 - 3) 


(6.3) 


t is in days and is zero at the vernal equinox. 4), is a mild time-lag introduced 

J 

in the troposphere to account for delayed ocean warmth. Trenberth used the fol- 
lowing values for his eight levels: 
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TABLE 6.1 


Z J 

P. 

J 

Zj (km) 

hj(day _1 ) 

d. 

J 

Tj 2 ( deg ) 

T*i (deg) 

Tj 3 (deg ) 

.511 

.60 

4 

0.080 

21 

-15.4 

6.54 

0 

1.204 

.30 

9 

0.059 

21 

-10.3 

5.97 

0 

2.120 

.12 

15 

0.019 

0 

- 2.24 

8.96 

1.70 

3.219 

.04 

22 

0.026 

0 

- 4.64 

11.80 

3.32 

4.605 

.01 

30.5 

0.050 

0 

- 5.39 

14.00 

4.57 

6.214 

.002 

43 

0.190 

0 

- 5.37 

11 .70 

5.45 

7.601 

.0005 

54 

0.240 

0 

- 1.78 

11.20 

4.89 

9.210 

.0001 

63 

0.200 

0 

0 

2.77 

1 .21 

The values of 

h. used in 

our model 

(whi ch 

are given 

in the next 

few pages) 


are similar to those of irenberth, in the troposphere particul arly . The values 
★ 

of T , however, have been rederived using Newell et al . [1972: The energy balance 
of the global atmosphere, in The global circulation of the atmosphere, pp. 42-90, 
Royal Meteorological Society, London 257 pp.; and 1974: The general circulation 
of the tropical atmosphere, Vol . 2, The MIT Press, Cambridge, MA 371 pp.] values 
of atmospheric heating rates. Thus 

T* ({!-) + V (6.4) 

P 

where q/Cp is obtained from Newell et al . and Tq represents the observed temper- 
ature distribution. 

We use 
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(6.5) 


T/u.t) - T* 2 P°(p) * TjVJ(u) + [TjV?(m) + T* 3 P|( U )] • single (t - $ .)j 

where T*. -1.1 (a®, + b° , ) . 

<j>j = 30 days, 

and there are twelve 30 day months in the model. 

The values of a*?, and b9 . are given in Table 6.2. The resulting values of 

J 

* 0 
T are similar to those used by Trenberth except for the inclusion of the P^(u) 

component. 


TABLE 6.2: Zonal heating parameters 


j 

Z 

■? 

(degrees ) 

0 0 
a 2 a 3 

•2 

■>? 

(degrees/day ) 

.0 ,0 

b 2 . b 3 

h° 

b 4 

32 

0 

(5.5 

-11.5 

0.5 

-1.5 

.20 

.06 

.08 

-.24) 

31 

.405 

4.5 

-10.0 

0.1 

-0.8 

.34 

-.20 

.15 

-.25 

30 

.811 

4.0 

- 8.9 

0.1 

-0.3 

.24 

-.53 

-.13 

.28 

29 

1.216 

3.7 

- 7.5 

0.2 

0.1 

.11 

-.15 

-.10 

.29 

28 

1.622 

2.9 

0.0 

1.3 

-0.1 

.07 

-.11 

-.05 

.13 

27 

2.027 

2.6 

7.4 

3.0 

-2.3 

.10 

-.10 

.00 

.02 

26 

2.433 

3.2 

8.4 

3.2 

-3.0 

.14 

-.20 

.02 

-.02 

25 

2.838 

4.4 

4.6 

3.3 

-2.1 

.17 

i 

-.20 

.03 

-.02 

24 

3.244 

5.2 

1.5 

3.2 

-1.4 

.19 

-.21 

.04 

.02 

23 

3.649 

5.6 

- 0.1 

2.7 

-0.8 

.21 

-.29 

.02 

.06 

22 

4.055 

5.9 

- 0.6 

2.0 

-0.6 

.29 

-.35 

.00 

.06 

21 

4.460 

6.0 

- 1.4 

2.0 

-0.5 

, 6 

-.35 

-.03 

.02 
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In addition to latitudinal forcing, the lower troposphere in the model is 

forced non-zonally. This forcing is applied at the three lowest levels only 

(levels 29, 30 and 31) and is based on Katayama's calculations of tropospheric 

heating rates. The heating is applied at each of the levels in the ratio P.h.. 

J J 

Thus, most of this forcing occurs at 667 mb. Fitting Katayama's computations 
1 2 

by P n (y) ancl P m (d) terms where n = 1, 2, 3, 4 and m = 2, 3, 4, 5, the heating 
may be divided into those terms which are even about the equator (the annual 
terms) and those which are odd (the seasonal terms). The even terms are found 
to be produced primarily by equatorial forcing whereas those which are odd are 
due to mid-latitude forcing. Since our model is quasi-geostrophic, it may not 
respond realistically to equatorial forcing. Although the even terms were in- 
cluded in the 6 wave version of the model (but excluded from Trenberth's model), 
they have been excluded from the 18 wave model. Early tests of this model pro- 
duced excessive tropospheric kinetic energy; the parameter changes made to pro- 
duce a more realistic simulation included the elimination of these "annual" terms. 
The seasonal terms are given by 

( 6 . 6 ) 


* 

T j 


P A 


sin lin (t 


31 360 

l p j h j 

j=29 J J 


l l (tc^cos l\ + is^sinZA) P^(y) 
J n=£+l , 1 = 1 n n n 

1+2 


where X is longitude. 


Values of tc 


l 

n 


and ts 


l 

n 


are given in Table 6.3 
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TABLE 6.3: Non-zonal tropospheric forcing 


l 

n 

< 

TS ^ 

1 

2 

1.12 

8.05 

1 

4 

-2.09 

-1.48 

2 

3 

-8.57 

2.15 

2 

5 

-1.23 

2.14 


In the stratosphere, Trenberth's formula empirically represents all of the 
main types of radiation effects: 

a. Short wave absorption by ozone and oxygen (20 km and higher) 

b. 9.6 micron absorption and emission by ozone (20-30 km) 

c. Infra-red absorption and emission by CO^ (all heights) 

In our model we will explicitly compute radiation of type a, as described in 
Chapter 5. This means that Trenberth's formulation must be changed for Z higher 
than 20 km (Z > 2.9). Type c can be represented by a simple Newtonian law 

f£= - h(p)T . (6.7) 

The first reported runs of the model (Cunnold et al . , JAS , 32 , 170-194, 1975) 
used values obtained by Dickinson (A method or parameterization for infra-red 
cooling between altitudes of 30 and 70 km, JGR , 78 , 4451-4457, 1973). He ob- 
tained values of h(p) by careful line integration of the CO and 0^ bands for a 
standard atmosphere, followed by a similar computation in which T was varied 
slightly. (Physically, this is satisfactory for the circumstance of cooling to 
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space - it does not include, however, the major effect of b above, which comes 
from absorption of radiation emitted by the ground in the water-vapor window 
region). Dickinson found values as follows: 

Z h(day _1 ) 


11.4 

.016 

10.0 

.062 

9.1 

.125 

8.7 

.172 

7.9 

.200 

7.4 

.220 

6.9 

.212 

6.1 

.135 

4.8 

.080 

3.9 

.060 


The values of h currently being used are given as h-j in Table 6.4. Between 
Z = 4.8 and 8.7, they match Dickinson's values. They, however, exceed Dickin- 
son’s values above Z > 8.7 and are smaller at Z = 3.9. The higher values at 
high altitudes are supported by the work of Kuhn and London ( J . Atmos . Sci . , 

1969) while the values at low altitude are consistent with the values of h given 
in Table 6.1. Arguments in support of the latter values are given in Phillips 
(Tropospheric and lower stratospheric heating, internal memo, June 28, 1973). 

The support is provided by the radiation calculations of Manabe and Strickler 
( J . Atmos . Sci . , 21 , 364, 1964) and by experience with numerical models of the 
troposphere (e.g., Bushby and Whitelam, QJRMS , 87 , 374-392, 1961). 
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TABLE 6. 
model 

.4: Two profiles of Newtonian 

calculations . 

cooling efficient used 

in the 

Model 

level 

Z 

h i 

(day -1 ) 

h 2 

(day x ) 

1 

12.6 

.06 

.06 

2 

12.2 

.03 

.08 

3 

11.8 

.09 

.09 

4 

11.4 

.10 

.10 

5 

10.9 

.11 

.11 

6 

10.5 

.13 

.13 

7 

10.1 

.14 

.14 

8 

9.7 

.15 

.15 

9 

9.3 

.16 

.16 

10 

8.9 

.17 

.17 

11 

3.5 

.19 

.19 

12 

8.1 

.20 

.20 

13 

7.7 

.21 

.21 

14 

7.3 

.21 

.21 

15 

6.9 

.19 

.19 

16 

6.5 

.17 

.16 

17 

6.1 

.14 

.14 

18 

5.7 

.12 

.12 

19 

5.3 

.11 

.10 

20 

4.9 

.09 

.09 

21 

4.5 

.07 

.08 

22 

4.1 

.06 

.07 

23 

3.6 

.04 

.06 

24 

3.2 

.03 

.05 

25 

2.8 

.02 

.05 

26 

2.4 

.02 

.05 

27 

2.0 

.02 

.05 

28 

1.6 

.03 

.05 

29 

1.2 

.05 

.06 

30 

0.8 

.07 

.07 

31 

0.4 

.10 

.09 

32 

0.0 

.14 

.14 


6-7 


The formula used in the computations to generate the values of h-j given in 
Table 6.4 is 

h = 0.22 - 0.03 (Z - 7.4) Z > 7.4 

3/2 

= 0.02 + 0.20 { Z 7.4 > Z > 2.5 

2 

= 0.02 + 0.12 i 2 ' 5 z l - } 2.5 > Z > 0 . 

Heating due to the 9.6 micron 0 3 band has been computed by Dopplick (Ph.D. 
thesis, M.I.T., 1970) using typical observed values of 0^. It gives a heating 
rate more or less independent of season with a peak value of 1.1 deg day ^ at 
the poles, centered at an altitude of about 25 km (Z = 3.5). This may be repre- 
sented most simply by 

= . M P 0 ( u ) [i . ( Z - 3.5) 2 ] deg/day (6.8) 

for 2.5 < Z < 4.5 and zero for Z values outside this range. As such, this heat- 
ing does not depend on 0 3 or T. This is not too critical, however, since the 
physical dependence on T is primarily due to ocean surface and cloud top temper- 
atures, and any moderate variations in 0^ from that used by Dopplick will mostly 
shift somewhat the height Z = 3.5 vertically (i.e., the existing 0^ is ample to 
absorb all this upwelling radiation). 

We therefore use the following strategy for heating computations. Define 

Method I. 

a. Explicit computation of short wave heating 
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b. Infra-red cooling with the Newtonian cooling formula (6.7 

c. Formula (6.8) for window radiation for 2.5 < Z < 4.5 

Method II. 

Trenberth's formulas 


These methods are combined as follows for different height ranges. (Z 
4.5 corresponds roughly to heights of 17.5 and 30.3 km). 


A. Z > 4.5: Only Method I. 

B. 2.5 < Z < 4.5: Weighted average of I and II as follows 


f Z - 2.5) 


x Method I + 


(4.5 - Z] 


x Method II. 


C. z < 2.5: Only Method II. 


2.5 and 
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7. Spectral form of the equations 

We define spectral solutions at arbitrary level j in the form 


*j = £Vj Y a (X - ll) 

Cj = Z5 a>j Y a (X,u) 

“j " 

T. = Vi , -Y (X,y) 
J L a J or 
a 

q j = ^a'J Y a (X ' u) 


(7.1) 


and for the trace substance equation 


*j = Ix a .jY a (X,u) 

a 

G . = IG , . Y (X,p) 
j L a J a 


(7.2) 


In terms of longitude (X) and latitude (p), we have defined members of the 
complete set of orthogonal spherical harmonics in (7.1) and (7.2) using 


M X 

Y (X,y) = e a P Jm) 


a 


a 


(7.3) 


wi th 


a 


n + il 
a a 


(7.4) 


denoting a vector index of planetary wave number t and degree n^. The P a (p) 
are Legendre polynomials of rank and degree given by a. Normalization of the 


7-1 


spherical harmonics is such that integration over the unit spherical surface 
(s) yields the orthogonal property 


★ 

Y Y q 

. a 3 


ds = 4tt6 


a’B 


(7.5) 


Complex conjugate values are denoted by an asterisk. Another useful property 
of the set of spherical harmonics is that they satisfy the differential equation 


V 2 Y = -c Y 


a 


a a' 


c = 
a 


n (n +1 ) 
or a ‘ 


(7.6) 


The complete set of orthonormal Legendre polynomials as used in (7.3) are de- 
fined such that 


P * = P 
a a 


(7.7) 


and all P have been normalized such that 
a 


+1 
- -1 


P P Q 
a 3 


dy = 26 


a, 3 


(7.8) 


We now want to substitute solutions (7.1) and (7.2) into the non-dimensional 

forms of our model equations, multiply through with a member of the orthogonal 
★ 

set (say, Y^), and integrate the resulting relationships over the unit sphere. 
Application of this procedure to the vorticity equation (3.9), for example, 
yields the desired spectral form of this equation, 
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,j 


, j ~ \,J + c Y _ £ [fr-l^Y-e.j+l - 
-(^Vej - A. [<7?T>«Y + e.j + 1 - 


Y+e 


“ ( FT )W Y+£,j + ( ??T )F Y,j+l ‘ ( r-l )F y,j + 


(2" c r ) A H) j 


”Y >J 


(7.9) 


in which, over the unit spherical surface s, 


dC 


Y.J - 


dt 


i l ip ■ 
Y Y.J 




1_ 

4ir 

1_ 

4ir 

1_ 

4tt 


34 


i Y* ds 


at y 


0(* f uJYJds > L 


J . Y 


diP, 


ax y 


■Y*ds 


J (ip . r. )Y*ds (See Appendix A) 
J) J Y 


(7.10) 


_X_ 


w 


y 


Ve Y-e.J 


Ve Y+£ ’ J 


1_ 

4tt 


[V-uVL(W.)]Y*ds (See 

J Y Appendix B) 


P = — 

Y.J 4ir 


F.Y*ds 
J Y 


(2_C Y )A H)jS,J ~ 4ir 


A H)j [V> j +2V 2 tj; j ]Y*ds 


Similarly, the thermodynamic energy equation (3.19), the trace substance equa- 
tion (3.21), and the thermal wind relationship (3.18) reduce to the spectral 
forms 
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dT 


dt ' ” B y ,j " S j W Y,J " A H)j C y T y,j + C Cp8n J a z 3 q Y,j 


id = _ B ix) . (^)w Y jJ . + ( r ^)G Y>j - C r li ) G Y ,j 


" A H ) i C vX v • + 


1 

n_ (dji) ' 

n m K dt’c 

e 2n 

J s 

m 


Y*ds 
; y 


H)j Y1 >j 4tt J 

VZ c T . = -D (tj; . , -<J> . ) + E(iJ; , . , , • ) 

Y Y »J Y Y-e.J-1 Y-£»J Y+e>J-l Y+£>J 


(7.11) 


where, for example, 


dT . 

Y,J 

dt 


1_ 

4tt 


9T i 

— ^Y*ds 

.at Y 


dt 


»J _ 


1_ 

4tt 


d Xn 


.dt y 


•Y*ds 


c y = } 

Y Y > J 4ttJ 


(-V 2 T.)Y*ds 
J Y 


B Y ,j = 8t? > T j) Y * ds ( See Appendix A) 


b y j = ,Xj)Y*ds (See Appendix A) 


D ii> . - E ii) , . = - — | 

Y Y-e>J Y Y+^.J 4 ttJ 


[V*yVt|j .]Y*ds (See Appendix B) 


r y Y 



In addition, we want to determine the spectral form of (1.6) relating the verti- 



! 


cal component of relative vorticity (c) and the streamfunction (^). It can be 
shown that 


’Y.J 




(7.13) 


or 




(7.14) 


provided that in (6.14) we stipulate y/O+iO (i.e., c^O). 

The spectral relationships (7.9), (7.11), and (7.13) [or (7.14)] along with 
definitions (7.10) and (7.12) form a complete set of equations for solution. 
However, it is not convenient to attempt to integrate the model in this form as 
there is no explicit relationship determining the vertical velocity field repre- 
sented by W. In order to define W, we want to alter the thermal wind relation- 
ship in (7.11). This development is contained in the next section. Furthermore, 
specification of the truncation limits to be used for series solutions (7.1) and 
(7.2) have not yet been establ ished' and will be discussed in a later section. 
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8. Determination of W in the dynamic equations 


In order to obtain an explicit description of the vertical motion fields in 
our model atmosphere, we insert (7.14) into the thermal wind equation of (7.11) 
and differentiate w.r.t. time to get 


AZ c 


dT. 
Y dv 


YvJ 


Y 


Y-e 


f d S-e,j-l _ d S-e ,j 


dt 


dt 


E y _ d S + £ > J 

: , dt dt 

Y+e v 


(8.1) 


for all levels j = 2,3, J-l. We note that (8.1) does not apply for the 

cases y = O+iO. Furthermore, for notational purposes, we will stipulate that in 
(8.1) and all future relationships, terms which require y~£ = 0 + i0 or n _<£ 
do not exist. This applies equally to cases in which y+z. is not contained within 
the specified model truncation limits. 

Let us now define 


a Y,j E ' Wj-1 “ Vd-1 + \,J ’ 

+ ^FT^ F y,j ~ ^FT^yJ+I + 

A H)j ( 2 " C Y ) ^Y > j-l ’ C Y>j) 
b Y,j E " B Y,j " A H ) j C Y T yJ + [ rW ]q Y,j 
such that using (7.9) we can write 
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d S,j 


i D 
1 Y 


^ [ V E ,i-l - (r+ 1 ) Ve,j + 


+ r Vc,j + l ] + I?^TTE± [W 


Y+e.j-l 


" (r+ 1 > Ve ,0 + rt' Y+e>j+1 ] 


and, the thermodynamic energy equation of (7.11) reduces to 


dT , 

Y»J = h - s w 

dt b y,j J Y > j * 


Inserting solutions (8.3) and (8.4) into (8.1) has the effect of eliminating 
the time dependence of (8.1) and at any given time we have 


AZ C Y b Y,J " AZ C Y S j W Y»J C Y _ £ S Y-£,j “ C y+£ VeJ •- 

" ('"-I) c^c— ^ W Y-2e,j-l " ^ r+1 ^ W Y-2e,j + rW Y-2e,j+l^ + 


i |E 0 EDY 
+ - T - -X-£ Y + - O+e 


T^ry Jc^C^ + C^C^j CW Y,j-l ~ (r+1 > W Y,j + rW Y,j + l^ 
T-l) c y +£ c y +2 £; ^Y + 2c j j-1 ” ^ r+1 ^ W Y+2c,j + rW Y+2e,j+l-^ 


or, if we define 
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n . .= (r-1 ) 
Y 5 J V ' 


JL 


C Y-e C y ^ £ c y c y+£ 


a . - AZ b 
Y + e,J Y >J 


m \ D D 
:( 1 ) = Y-e Y 


Y 


'Y-2e C Y-£ C Y 


:(2) - 1 


Y 


fE D ED. 
Y-e Y + Y Y+e 


'Y 


'Y-e Y+e 


(8.5) 


* ( 3 ) _ E Y E Y+e 


c.,c . c 


Y y+e y + 2 


= (r " 1 )AZ s j 


the W-equation can be compacted to 


r f O) w + f (2)u + f^2)u i _ 

^ Y Y-2e,j-l Y t » j-1 Y W Y+2e,j-l-* 


- < r+ ’> + f y\,i + f i 3, Vz.J ] + 


+r[f^^W ? . , + f^W ,j+l + f| 3 ^ 

Y Y~2e,J+l Y Y Y 


Y+2e,j+l 


- 


( 8 . 6 ) 


- a-W . = y . 
J Y.J Y i J 


in which from (1.17) we represent the boundary conditions as 


W , = 0 
Y > 1 

W , = H 

Y,J Y 


(8.7) 


H =1-| 
Y 4tt J 


. j( Vi, h/H o) Y ; ds 
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To prepare (8.6) for inversion we want to take note of certain properties 
of the equations in order to reduce the calculation to a finite set of simple 
matrix solutions. Inspection of (8.6) shows that the equations uncouple accord- 
ing to planetary wave numbers, Z . In addition, within each planetary wave the 
equations contain two independent sets; one of even vector elements (n + Z^ all 
even) and the others of odd vector elements (n + Z all odd). Thus, to facili- 
tate ease of notation, let us define some new sets of indices to be applied to 
(8.6) by first denoting a maximum planetary wave number, L, for a given spectral 
truncation as 


L 


^max 


( 8 . 8 ) 


so that we can designate K independent sets of matrix equations using index k 
where 


k = 1, 2, 3, .... K; K = 2(L+1). 


(8.9) 


For a given matrix set we will determine k by designating 

2 Z + 1 for even vector sets 

Y 

2 (Z +1) for odd vector sets 
Y 




( 8 . 10 ) 


Furthermore, within each of the K matrix equation sets it is useful to designate 
an element index, b^, where 


b^ - 1 , 2 , 3 , . . . , 


( 8 . 11 ) 


Thus, for a given matrix set designated by the subscript k we devise the b^ 
indices as follows: 
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(1) For k odd (even vectors) let 


^k ” n k^max 


( 8 . 12 ) 


for which we consider only n^ from the set n^ + even. Then the value for an 
individual is determined from 


(8.13) 


where we ignore values of outside the range indicated in (8.11); i.e., when 
k = 1, n-j = 0, and £-| = 0 we do not include the value b-j = 0 which designates 
the nonallowable equation of (8.1) in which y = O+iO [see comments following 
( 8 . 1 )]. 

Similarily, 

(2) For k even (odd vectors) let 



N k = 


n k^max 


(8.14) 


in which here we consider only n^ from the set n^ + odd. Then, we have 


b k = 


n k - *k + 1 


B k = 


N k - ^k + 1 



(8.15) 


At this point we want to note an additional property inherent in the spec- 
tral W-equations represented by (8.6). That is, from definitions contained in 
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(8.5) and Appendix 8 we can show that for any given k, 


E E . 

f (3) „ Y k Y i< E 
D. c c c 

k Y k Y iE 2e 

(8.16) 

= V V £ f (l) 

CY k Cy t l: Cy t 2e 


We are now prepared to convert (8.6) to matrix form. To do this we first 
define tri diagonal matrices V ^ as 


V 


k 


fj 2 ^ fj 3 ) 0 


0 


f ( 3) f (2) f (3) 

r l r 2 T 2 


6 





0 


f (3) >(2) 

V 1 B k 


(8.17) 


where we have made use of (8.8) - (8.16). We note from (8.17) that not only is 
each tri diagonal but it is also symmetric. In addition, it can be shown that 
every principle minor determinate of is posi ti ve and thus V ^ can be said to 
be positive definite . These properties will be discussed in more detail below. 

To complete the conversion of (8.6) to matrix form we define vectors 
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w l,j 


t \ 


W 2,j 


, J 

w k,j = 

V 

’ R k,j = 



V 

V. j 

k 

V. ■/ 


such that (8.6) can be written in the matrix form 


Vk.J-1 ■ + rt ’k w k,o>1 - a j w k,j ■ R k,j ; 

j = 2, 3, 4, , J-l for each k = 1, 2, 3, , K 


(8.19) 


We wish to modify (8.19) through di agonal izati on of each V^. However, 
since each tri diagonal V ^ is real , symmetric and positive definite , we know that 
all eigenvalues of Pj are real and positive . Also’, the sets of eigenvectors 
associated with these eigenvalues are orthonormal . Thus, if V. is an MxM matrix, 
there exists a set of real positive eigenvalues (A, ) with p = 1, 2, 3, ..., M 

K p 

associated with V. and M sets of orthonormal eigenvectors q with s = 1, 2, 3, 

K r 

. . . , M. If we let represent the matrix of eigenvectors associated with the 
the set (A k ) and matrix V^, we have 


q n 

q 12 

... q ls 

••• q lm 

q 21 

q 22 

••• q 2s 

••• q 2m 

q pi 

q p2 

- q ps 

q pm 

q ml 

q m2 

. . . q 

. . . q 
M mm 


8-7 


such that 


Q k Q k ~ Q k Q k = 1 


( 8 . 21 ) 


where I is the unit matrix and ( ) denotes transposition. Define 


f(x k ), 0. 


A k " 


0 -U k ) 2 




(XJ 


m 


( 8 . 22 ) 


were then we know 


and 


P k Q k Q k A k 


Vk Q k ~ Q k Q k A k " A k 



We now want to expand the vector W. . in (8.19) in the form 

K 5 J 


(8.23) 


w k,j - Vk,j ; v kj = Q k w kj 


(8.24) 


where we note that V. ■ is also a vector. 

K , J 


Inserting solutions (8.24) into (8.19) and multiplying through with gives 

WkYj-1 - < r+1 >Wk v k,j + r \Wk,j*i - 

- “jWk.j * V*k,j 


or, from (8.23), we can write 
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(8.25) 


A k V k,j-l ' ^ r+1) A k + a j^ V k,j + rA k V k,j+l " Q k R k,j 


Now, we know that there exists an inverse 


A, 


i/(x k ) 1 0 

o l/(A k ) 2 


i/(x k ) t 


l/(xj 


k'm 


such that 

\\ = 1 • 

Thus, if we multiply (8.25) through with A”^ , (8.25) reduces to the form 

V k,0-1 - C(r+1)I + °A 1] \,j + rV k,j+l = A k\ R k,j 
where for each k = 1 , 2, 3, . . . , K we have j = 2, 3 

\,j - ' [(r+1)I + 0 o A k 1] 
and 



: 2, 3, 4, . . . , j-1 . 

We now 

A'] 


k.i (for J “ 2) 



(for 3 < j < J - 2 ) 

l k' Q k R k,j-i - rV k,j (for j ■ J - 1) 
Using (8.29), (8.28) transforms to the set 


(8.26) 

(8.27) 

(8.28) 
let 

(8.29) 
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S k,2 V k,2 + r V k,3 ~ R k,2 ^ for J ' 


V k J-l + 5 k ,j V k,j + r V k,j + 1 = R k,j {3 ^ J ~ 2) ^ 8 - 3 °) 


V k,J-2 + S k,J-l V k,J-l R k,J-l ^ for j J_1 ) 


in which from (8.24) and the boundary conditions of (8.7) we see that in (8.29) 



Vi - 0 


(8.31) 


V k,J - «k W k.J = «k H k 


We see that for each k the system (8.30) is tridiagonal in j and thus submits 
readily to solution provided certain provisions are met (see Appendix C for 
details). Briefly, to invert (8.30) we first define 


and then let 


u k,2 

- S k,2 

(for j = 2) 

) 

U k,j 

E (S k,j “ r U k,j-l) _1 

(for 3<j<J-l) ^ 

> (8.32) 

V k,j 

E “ r u k,j 

(2<j<J-l) ^ 

1 

y k,2 

= u k,2 R k,2 

(for j = 2) ~~J 



Yu. = u. , (R. , - y u , ,) (for 3<j<J-l ) 


k,j k, j v k,j ^ k , j - 1 

Solutions to (8.30) thus appear as 


. (8.33) 


8 -'. 0 



(for j = J-l ) 


(8.34) 


V k,J-l = y k, J-l 

V k,j = V k,j V k,j+l + y k,j (for j = J-2, J-3, . .., 2) 

provided all u R . in (8.23) exist and are finite. Vectors W R ^ 
tained from (8.24). 



are then ob- 
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Appendix A . Evaluation of the nonlinear Jacobian terms using transformed 
fields . 

On the unit sphere, the Jacobian of arbitrary horizontal global scalers 
A and B is given by 


i / « D \ _ 3A 9B 3A 3B 
J ^ A,Bj 3A 3y ' 3y 3A 


(A. 1 ) 


where X is longitude and y is the sine of latitude. Expanding A and B in terms 
of spherical harmonics, we have 


A = Ya Y (A,y) 

L a a ' 
a 

8 = ^ b a Y a (X ’ y) (A ' 2) 

a 


in which the special properties of the orthonomal spherical function Y a (A,y) are 
outlined in (7.3) - (7.8). 

Rather than evaluate the Jacobian terms using the interaction coefficient 
methodology applied in the past for our lower order models, we will transform 
the terms in (A.l) to a special grid system devised to permit exact Gaussian 
guadrature integration of the nonlinear Jacobian terms. Details of this numer- 
ical procedures to be used in this process are contained in Eliasen et al . , 
(1970). Briefly, however, we can assert that the selected grid system must be 
located on latitudes which form the "zeros" of the "highest order plus 1" Legen- 
dre polynomial contained in the Jacobian term for the particular truncation 
chosen. In the case of triangular spectral truncations (which we are currently 
using), we know that the number of the "Gaussian" latitudes needed, K is 
given by 


A-l 


K max 1 ' 5 Anax + 1 ’ 


(A. 3) 


where l is the maximum longitudinal wave number allowed. Thus, for a trunca- 

ITla X 

tion triangle which contains planetary waves, 0 - 18, the formula requires that 
28 latitudes be used to represent the minimum horizontal grid and that the 1 a t i - 

O Q 

tudes are to be located on the zeros of the Legendre polynomial Pq . A similar 
procedure is used to specify the longitudinal grid points. 

To evaluate (A.l) using the transform technique, we first define 


= a Y (A ,y) 

9a l a a a' 
a 


3B 


3u 



iv dl >) 

dy 


3 9A 
3u 



U AdP (y) 
a a VH/ 


dy 


(A. 4) 


Q = 


3B 

9A 


a 


where then (A.l) becomes 


J(A,B) = MxN - PxQ . (A. 5) 

The quantities M, N, P, Q are evaluated on the Gaussian grid by performing the 
sums indicated in (A. 4) for each grid point. Details of the numerical methodol- 
ogy used for these transforms are given in Appendix E. The indicated products 
and differences indicated in (A. 5) are performed at each grid point and the 
resulting Jacobian field is transformed back to the spectral domain and trun- 
cated to conform with the resolutions of the A and 3 spectrums. 
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Reference 


Eliasen, E. , B. Mackenhauer and E. Rasmussen, 1970: On a numerical method for 

integration of the hydrodynamical equations with a spectral representation 
of the horizontal fields. Report No. 2, Institut for Teoretisk Meteorologi, 
Kobenhavns Universitet, 37 pp. 
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Appendix B . Spectral representation of divergence terms of the general form 
V • yVA . 

In terms of spherical operators on the unit sphere in which X is longitude 
and p is the sine of latitude we have 

V • pVA = Vy • VA + yV 2 A 


(B.l ) 


- o-u 2 ) + uv 2 a 

op 

in which A is an arbitrary horizontal global scalar expandable in the form 

A = la Y (X,u) . (B.2) 

a 

Properties of the orthonormal spherical functions Y (\,u) are outlined in (7.3) - 

ot 

(7.8). Insertion of solutions (B.2) into (B.l) yields 


V • yVA = (l-y 2 )^ 1 ^ 


a 


= 7a e 

L a 

a 


i Iq\ 


dy a 

dP 

‘ pC a P a 


c = n (n +1 ) 
a a' a 



(B. 3) 


But, if we define 


N = 
a 


(2n +1 ) (n -i ) 
v a a a 


tV2 


(n +1 ) ! 
v a a 


(B.4) 


c = 1 + id 
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then we know from the Legendre differential and recurrence relations (for exam- 
ple, see Jahnke and Emde, 1945) that 


dP N 

0-d 2 )t^ = - n ~U p ~ + ( n „ + ON P 


dy 


a" a 'a a N a _ e a-e 


and 


p . Wl\ p + ( WA_p 

u a “T2T+TT N a+£ a+e {ZT+TJ N a _ £ a-e 



(B.5) 


Then, using (B.5), we can show that 


dP 

n 2 \ a D - x ' "a' '"a a' a p 
( u ^dy uC a a (2n^+l ) a-e 


(l-n 2 )(n +1 ) N 


a 


a-e 


n (n +2) (n -l +1 ) N 
a a a a a p Cpr'i 

C2F+T5 Ca+e 1 (B - 6) 


a 


a+e 


We now insert (B.6) into (B.3), multiply through using Y* / 4 tt , and integrate 
over the unit sphere to get 


T_ 

4ir 


2 it f 1 

(V-yVA)Y*dydX 
-1 Y 


I a 


a 


CL 

l =1 
a y 


(1-n 2 ) (n +1 ) N 
' a a a a 


(2n +1) N 
a a-e- 


r 1 


1 

2 


P P d 
_1 a-e y y 


I 


a 

l =1 
a y 


n (n +2) (n -l +1) N 
a a a a a 


(2n +1) 

a 


a+1 


1 

2 


-1 


P + „Pydy 
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(1 -n 2 ) 
Y 


(n X )(n -l ) 
Y + Y Y y J 


(2n Y -l)(2ry-]) 


1/2 


Y-e 


n Y (n y +2) 


( Wl)(n Y -yi) 


(2n Y +l)(2n Y +3) | “y+e 


D 3 -Ed 

Y Y-e Y Y+s 


(B. 7) 


where we have defined 


° Y = 0-"*) 


( n Y + ^Y^ n Y~^Y^ 
(2n -l)(2n +1) 


n 1/2 


E y = n Y ( n y + 2) 


Y 


(VV 1)(n YV ]) 

" (2n'+l)(2n +3) 


(B.8) 


Y 


1/2 


A special case of (B.7) occurs when we consider scalar B in which 


B = V 2 A 


(B.9) 


where similar to (B.2) we can expand B in the form 


= Zb 

a 


Y (A,u). 
a a 


Then, from (7.6), we know that 


b = -c a 

a a a 


and, in terms of coefficients b , (B.7) becomes 

a 


(B.10) 


(B.ll) 
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(B. 12) 


T_ 

4tt 


■2ir r 1 

(V*yVA) Y*dydX = 

J o j -i Y 


= V b 
Y Y-e 


D 

Y 

c 

Y-e 


b 

Y-e 



E b , 

Y Y+e 


in which we have defined 


V 

Y 



Y-e 


E 


Y 



(B- 1 3) 


provided that in (B.12) we ignore terms in which c = 0 (i.e., n = 0). 

i Y 

Further, for both (B.7) and (B.13) we must stipulate that all terms calling for 

any a , a , , b or b , outside the range of the particular spectral trun- 
y-e y+e y-e y+e 

cation chosen must also be ignored. 


Reference 

Jahnke, E. and F. Emde, 1945: Tables of Functions . Dover, New York, 306 pp. 

plus tables. 


B-4 



Appendix C . Solution of a tridiagonal set of equations . 


Suppose we have an equation set of the form 


aX -,+bX+c X,, = R 
Y Y-l Y Y Y Y+l Y 


Y = 1, 2, 3, .... r 


(C.l) 


where we must have 


a-| = 0 

c = 0 
Y 

That is, in matrix form we can write (C.l) as 


AX = R 


with A being tridiagonal of the form 


A = 



C 1 


0 


a 2 


0 


b 


2 


C 2 



(C.2) 


(C.3) 


(C.4) 


For solutions we define 


C 1 = l/b 1 

C = 1 / (b -a c ,C 2<y<r 
Y Y Y Y“ ' Y” ' - - 


(C.5) 


C-l 



B = C (R -a B ,); 2<y<T 
Y Y Y Y Y-l 


(C . 6 ) 


Then, the solutions appear as 

x r = B r 

X = D X , -,+B ; Y = r-1, r-2, 1 

Y Y Y + I Y 

provided all C in (C.5) are finite. That is, if 

b 1 t 0 

b t a c -,C n 
Y Y Y-l Y-l 



(C. 7) 


(C.8) 



Appendix D . Computation of the weight functions for Gaussian quadrature . 


/ 

We consider the set of complete orthogonal Legendre polynomials, P (p), in 
which l = 0, ±1 , ±2, ... and n = 0, 1, 2, ... . We define this set, according 
to (7.8), to be normalized such that 

,1 

P>)^(u)du = 2 Vn . (D.l) 

-1 


where p is the sine of latitude or equivalently, the cosine of colatitude, <j>. 

Now in order to expand an arbitrary function of latitude, say f(p), in terms of 
the set of Legendre polynomials we let 

f(u) = nfM(u) (D.2) 

In n n 


from which the coefficients, f^, are obtained through application of (D.l) such 
that 



2 


,1 

f(u)P^(u)dp • 
• -1 


(0.3) 


However, to be able to transform at will between spectral and grid point space, 
it is necessary to represent f(p) at a number of discrete points, p^, in which 
k = 1, 2, 3, ..., N with N being the total number of points lying within - 1 <p< 1 - 
Thus at each latitude point, (D.2) becomes 

f K>'^ P nK> ■ (°.4) 

In 

This means that in order to determine coefficients f^ we must evaluate the inte- 
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grals in (D.3) numerically and at the same time maintain the orthogonality prop- 
erties of the discrete polynomials representation in (D.4). For this purpose, 
integrating by quadratures, we introduce a set of Gaussian weight functions, w^, 
such that 




p£(u)P^(v)dy 


(0.5) 


-1 


and the numerical analog for (D.3) becomes 






(D.6) 


The remainder of this Appendix is devoted to the method of evaluation of the 
Gaussian weights, w^. 

Because we know that any given Legendre polynomial, p (p) , can be repre- 
sented by a finite series in p of at most degree n, we can expand 


or 


and thus, 


o o i 

- f 0 v 

P nK> P n'Kl " T 0 b iKf 


(D.7) 


A 

P^(p)P^(u)du 

-1 


n+rT 

l \ 

i=0 1 


A 

y^p . 

-1 


(0.8) 
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Integrating (D.8) by quadratures using (D.5), 


1 

-1 



N 

= Iw 
k=l 


n+n' 

k l b i 
K i=0 


[yj 1 


Equating (D.8) and (D.9) we have 


n+n" 

l b i 

i=0 1 


. N n+n" 
y^y = l w. I b,[yj 
k=l K i=0 1 K 


(D.9) 


(D.10) 


and thus for any i such that 0 < i < n + n" it must hold that 


rl 

y^y 

• -1 


l w kty k ^ 

k=l K K 


(D.n) 


We see from (D.ll) that if we choose the number of latitude points, N, such 
that N-l = n+n" then utilizing all i = 0, 1, 2, n+n we can form a set of 
N equations containing N unknown quantities, w^, for inversion. However, in 
terms of colatitude, <j> , we can show that any cosj<£ (j is an integer) can be 
expanded in the form 


and 


cos j 



r .j <V2)-' r ,2m 
cosj4>l, = a jCv- k 3 + J, a 2m [>V 


m=0 



(D.12) 
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Then, inserting (D.12) into (0.11), 


rl 


a . 


1 J-l 


COSifpdy - 


1 ( i / 2 ) - 1 
aT I 

i m=0 


2 m . 

U dy 


-1 


1 N , (i/2)-l N 

= l w k cosi<!) k " T % a 2 m ^ w k^ 

a i k=l K K a i m =0 k=l K K 


2m 


(0.13) 


or 


l w.cosi^t, = cosi(j)dy 
k=l K K -1 


cos i sin(j)d<j) 


'0 


0 for i odd 


—7 for i even 

i -1 

i = 0 , 1 , 2 , ..., n+n " 


(D.14) 


where we have made use of (D.ll) to eliminate the second term on each side of 
(0.13). Again, as for (D.ll), we see that if we take N-l = n+n", we can invert 
(D.14) to obtain the Gaussian weights. 

As an example, consider N=3 where we select 4>-j=30°, <f>2=90°, and <J>2=150°. 
Then, from (0.14) we can construct the set (using i = 0, 1, 2) 



(0.15) 
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with solutions 


W 1 = w 3 = 4/9 
w 2 = 10/9 


(D.16) 


We note that the solutions (D.16) are symmetric in about the equator. If we 
assume such symmetry a priori then all equations in (D.14) involving odd values 
of i become redundant and we can write (D.14) over the integration interval from 
<j> = 0 to 4> = tt/2 as 


N+l 


l w.cos2i<j). 
c=l K K 


tt/2 


0 


cos2i(j) si ncj)d4> 


1 


i = 0, 1, 2, 


W^T * 

n+n 


(D.17) 


; N-l = n+n' 


Again, using the example used above in which N=3, <t>-| = 30°, and <t >2 = 90°, we 
N+l N-l 

have -tj— = 2 and = 1 giving the set 



with solutions 



(D. 18) 
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Furthermore, if we want to obtain w^ s for the entire pole to pole integration, 
we need only make use of the symmetry property 


Vl-k = w k + \,,r/2 w k 

which gives for our example 



Solutions (D.20) are identical with those of (D.16). 


(D.19) 


(D.20) 
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Appendix E . Transform methedolegy . 

The transform routines take advantage of a number of factors related to the 
particular configuration of our model which allows for a considerable improve- 
ment in speed and efficiency. For example, since the physical fields that are 
being represented in our model are all real quantities (as opposed to complex 
values), we can make use of the inherent complex nature of Fourier transform 
theory to transform two of our fields simultaneously. Further savings in compu- 
tational resources are possible by taking advantage of the odd/even nature of 
the Legendre polynomial expansions with latitude about the equator. 

Consider the series representation 


4»UiV) = l P a (u) (E- 1 ) 

a 

in which we want to be able to transform freely between spatial data as repre- 
sented by (X ,u) in (E.l) and coefficient data as determined by the {^ a > for a 
given truncation. 

Now, in terms of even and odd series (with respect to the equator), we can 

write 


#U>u) 


l 

a=even 


if; e i£aX P (y) + 
y a or ' 


l 

a=odd 


e i£ ^Pju) 


a' 


(E.2) 


For simplicity we define 


E.(p) - l Pn“Cu) ; 


-v 


n„ n n 
a a ct 


■t Z 

= ^ if) a P a (y) ; a = odd 


; a = even 


(E.3) 


°^> ■ S a V "a 


where then E^(y) is even and 0^(y) is odd (with respect to the equator). Thus, 
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for y > 0, 


iJ'U.u) = l [^(u) + 

^U,-y) = l [E^(y) - 0^(y)] e 1 '^ 

We now let (for y > 0) 

A £ (y) = (1+i) E £ (y) + (1-i) 0 £ (y) 

B £ (y) = (1+i) E*(y) + (1-i) 0*(y) 

in which denotes complex conjugation. We then can write 
A(A,y) = i|>(A,y) + i^(X ,-u) 

= Kn + i)E/(u) + (l-i)O.(u)] e ilX 

l 

- lA.(u) e UX 
l C 

B(A,y) = ip(X ,y) - i^(A,-y) 

= [[(l-i)E o(y) + (1+i )0o(y)] e UA 
l ^ 




(E.4) 



(E.5) 


(E- 6) 



Using the definitions and relationships developed above, we now turn to a 
description of the computational routines developed for the spectral to grid 
and grid to spectral transformations required in the model. 


E-2 


I. Spectral to grid transformations 

a. Subroutine SPCGD2 (DSPEC, DATA) 

This routine transforms spectral coefficients (stored in "DSPEC") to 
grid values ("DATA") at the Gaussian latitudes. 

The routine calls subroutines "SPF0R2" and "FFTFR2". 

Storage locations for the input coefficient matrix "DSPEC" and the 
output grid point matrix "DATA" are illustrated below. 

1. Storage arrangement for "DSPEC". Assume a matrix of complex 
spectral coefficients, ij > , in which the vector index "a" is as 
defined, for example, as in Appendix A (A-2). The index "NLEV" 
represents the total number of model levels. Then, 



2. Storage arrangement for "DATA". Assume ip now represents elements 
in a matrix of real gridpoint values. Here "NLON" and "NGRID" 
are, for a given model level, the number of longitude points along 
a given latitude circle and the total number of gridpoints, respec- 
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tively. Note that data for the conjugate (i.e., Southern Hemi- 
sphere) latitudes appear in the last 32 columns of the matrix. 
Thus , 


(LEVEL 1) 

^(^2 jU-j ) 


NLON ^ ^ LON 5 ^“1 ^ 

NLON+1 ip(X 1 ,u 2 ) 

NLON+2 


^(X*j sU-j ) 

4* ( ^ 2 5 ^] ) 


^ ( X-| 5 U 2 ) 
4 1 ( 3 ^ 2 ) 


33 . . 

. . . 32 + NLEV 


(LEVEL 1) 

(LEVEL NLEV) 


^(X-j > -Ui ) 

. . ^(Xp-up 


4^2 ~y 1 ) 

ip ( X^ > “U"] ) 



^ X NL0N ,U f ^ X NL0N ,_U 1 ^ 


'l 1 ( > "V*2 ^ 

^(x 2 ,-u 2 ) 


■ • * (X NL0N’- U 1 ) 

. . ip(Xi ,-u 2 ) 

• • ^(^2 5 ~^2 ^ 


2 *NL0N ^^MLON’^2 ^ 
2*NL0N+1 4j(X-j ,y 3 ) 


^ X NL0N ,U 2^ ^^NL0N’“ U 2^ 
U;(X-[ ,ia 3 ) ^(X 15 -u 3 ) 


NLON EQ 


^ , ^NL0N’"' J 2^ 

^(X-| > "U 3 ) 


NGRID/2 ^(^ nlon »U eq ) • • - ^ X NL0N’ U EQ^ ^ X NLON ,_U E(p ’ ’ ‘ ^ X NL0N ’ _u E(p 


Subroutine SPF0R2 (DSPEC, DATA) 

This routine determines Fourier coefficients (returned in "DATA") for 
each latitude using the spectral coefficients in "DSPEC". 

1. "DSPEC" is the input coefficient field as in subroutine "SPCGD2" . 

2. The complex Fourier coefficient output is returned in "DATA" in a 



form which combines the conjugate latitude values. Thus, "DATA" 
will contain the A^(y) and B^(u) fields defined in (E.5) as re- 
required for the grid value expansions of (E.6). In the following 
matrix table, is the maximum planetary (or longitudinal) wave 

number, LR1 = L^^+l , and "NLAT" is the number of Gaussian lati- 
tudes. All other quantities have been defined previously. The 
storage arrangement for A^ and is as follows: 


(See following page) . 
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2*NL0N ReCB-j (u 2 )> 



NLEV s' 33 

(LEVEL NLEV) H (LEVEL 1) 


Re{A Q (u-|)} Im{Ag(ii-| ) } 

Re{A-| (jj-| ) } J? Im{A-| (u-j ) > 


Re(A. (vi,)> Im{A, (y-> )} 
l MAX 1 ' l MAX 


Re{B, Cut ) } IMB. (p-.)} 
L MAX 1 MAX 


Im{B 2 (ui ) } 
Im{B-| (u-j ) } 
Im{A 0 (u 2 )} 
Im{A n (y 9 ) } 


Im{B-j (u 2 )} 


32 + NLEV 
(LEVEL NLEV) | 

I 

WA 0 (y,)( \ 




Im{ A, (u 1 )} 4^ 
U MAX 1 ! 


Xm{ B . (u-,)} 

L MAX 1 


Im{B 2 (u 1 )} 
Im{ B-j (u-j )> 
Im{A Q (u 2 )} 
Im{ A-j (u 2 ) } 


Itn{ B-| (u 2 ) } ^ 


NLAT/2*NL0N Ra{B,(upn)> 


Im{B-, (lirn ) ) 


I m { B -j (u^q) } 4 ^ 
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c. Subroutine FFTFR2 (DATA) 

This routine converts the Fourier coefficients obtained in "DATA" 
from "SPCFR2" to grid values (also stored in "DATA"), as described 
for subroutine "SPCGD2", using fast Fourier transform (FFT) techniques. 
The transformation can be reversed by changing the sign on the flag 
variable, "ISIGN". It thus is used as the first step in transforming 
grid data to full spectral coefficient data. 


II. Grid to spectral transformations. 

a. Subroutine GDSPC2 (DSPEC, DATA). 

This routine reverses the procedure of "SPCGD2" and produces spectral 
output in "DSPEC" from gridpoint data contained in "DATA". The storage 
arrangements for these two fields are as described in "SPCGD2". The 
program begins by moving one latitude (and i.ts congugate latitude) of data 
from "DATA" to "DSPEC" on a temporary basis and then calling subroutine 
FFTFT2 (DSPEC) with argument "DSPEC" in order to compute the A^ and 
fields for that latitude. 

b. "GDSPC2" then continues by combining the A^(y) and B^(p) fields (again, 
by latitudes) to obtain fields of E^y) and Og(y) as defined in (E.3). 
These are given (for l > 0) by 


E^(y) = -J-[A^(y)+B^(u)] + {-[A^(u)-B^(y)] 
0g(y) = -^[A^(y)+B£(y)] - ^{A^(y ) -B^(y) ] 



(E.7) 


and, for 1=0, 

Eq(u) = j[Re.{A 0 (y)} + Im{A Q (y)}] 
0 Q (U) = ^[Re-(A 0 (y)} - Im{Ag(y ) }] 



(E.8) 


where Eg(y) and 0g(y) are both real quantities. Due to limitations in 
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storage, however, the E^(y) and 0^(y) values are stored both in the 
"DSPEC" and "DATA" matrices. That is, it was found to be convenient 
to place the and 0^ coefficients for the first three planetary 
waves {1=0, 1, and 2) in the region DSPEC ( 1,1 07 ) - DSPEC(64,190) . 

All of the other coefficients (£ > 3) are placed in. "DATA" by latitudes. 
The storage arrangements for both "DSPEC" and "DATA" follow. First, 
however, we want to note, for ease of computation, that we can make use 
of the definition, for example, 

A f = Im{A} + iRo{A} 

so that (E.7) can be written as 

E^y) = .25 [A £ (y)+B*(y)+(A*(y)-B £ (y)) f ] 

0 £ (y) = .25 [A £ (y)+B*(y)-(A*(y)-B^(y)) T ] 

Define, NLATHF = (NLAT+1 )/2 and J2 = N LATH F+N LATH F such that the 
storage for "DSPEC" is: 



(See following page). 



v m 


For "DATA" we let J2 = 2*NLATHF*(L MAX -3) where then: 
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